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SECTION  I 
INTRODUCTION 

A basic  problem  In  the  two  dimensional  flow  of  an  ideal  fluid 
is  the  determination  of  the  stationary  flow  around  a body  or  bodies. 

The  flow  field  is  completely  specified  by  either  the  velocity  poten- 
tial function  $(x,y)  or  the  stream  function  t(x,y).  If  either  is 
known  the  velocity  field  can  be  determined,  the  pressure  coefficient 
can  be  calculated  from  Bernoulli's  equation  and  the  forces  acting  on 
the  body  or  bodies  can  be  evaluated. 

Both  the  velocity  potential  function  $(x,y)  and  the  stream 
function  fCx.y)  satisfy  Laplace's  equation  at  all  points  exterior 
to  the  body  or  bodies.  The  normal  derivative  of  ♦(x.y)  vanishes  on 
the  surface  of  the  bodies  and  the  stream  function  TCx.y)  is  a con- 
stant on  the  surfaces.  The  gradients  of  both  functions  tend  to 
finite  limits  at  infinity.  These  properties  characterize  the  func- 
tions 1i(x,y)  and  T(x,y).  Thus  the  basic  problem  reduces  to  that  of 
determining  a solution  of  Laplace's  equation  in  the  region  exterior 
to  one  or  more  contours  and  satisfying  conditions  prescribed  in  ad- 
vance (boundary  conditions)  on  these  contours  and  at  infinity. 

Functions  which  satisfy  Laplace's  equation  are  called  harmonic 
functions.  In  potential  theory  the  functions  of  an  extensive  class 
of  harmonic  functions  are  expressed  as  either  the  logarithmic  potential 
of  a single  layer  of  density  k(s)  or  as  the  logarithmic  potential  of 
a double  layer  of  moment  (density)  m(s)  distributed  over  the  contours. 
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In  this  report  we  use  each  representation  (separately)  in  an  expres- 
sion for  the  desired  solution  of  Laplace's  equation.  The  solution 
for  the  basic  problem, is  obtained  in  principle  when  the  density 
function  is  detetTnined.  The  density  function  is  the  solution  of  an 
integral  equation. 

The  velocity  tangent  to  the  airfoil  may  be  considered  as  the 
end  product,  for  the  determination  of  the  pressure  coefficient  is 
routine  after  the  tangential  velocity  is  known.  Thus  in  SECTION  II 
•e  express,  in  turn,  the  velocity  potential  function  and  the 

stream  function  f(x,y)  as  a linear  combination  of  certain  simple 
harmonic  functions  and  the  logarithmic  potential  of  a single  and 
double  layer  respectively.  From  these  representations  of  the  vel- 
ocity potential  and  the  stream  function  we  determine  expressions  for 
the  velocity  tangent  to  the  contour.  We  see  there  that  the  tangential 
velocity  can  be  obtained  conveniently  in  two  different  ways. 

The  tangential  velocity  is  determined  as  the  solution  to  a 
Fredholm  Integral  equation  of  the  second  kind  or  depends  upon  the 
solution  to  such  an  equation.  In  SECTION  III  we  consider  the  re- 
duction of  these  integral  equations  to  systems  of  linear  equations 
and  give  the  reasons  for  our  preference  of  method  for  achieving  this 
reduction.  In  SECTION  IV  we  describe  and  discuss  our  numerical  e'-- 
periments  and  give  our  conclusions. 

The  methods  described  in  this  report  are  Integral  equation  meth- 
ods. Integral  equation  methods  have  been  carried  out  successfully 
In  past.  Indeed,  practically  all  methods,  even  conformal  mapping 


methods,  have  an  Integral  equation  in  their  development.  There  are, 
however,  some  variants  in  the  present  report  which  we  believe  are 
worth  notice. 

When  the  contour  is  a circle  the  kernel  in  the  integral  equation 
reduces  to  a constant.  In  this  case  the  correct  solution  of  the  in- 
homogeneous problem  is  obtained  in  one  Iteration  step,  the  reason 
being  that  here  all  eigenvalues  of  the  kernel  are  zero  except  for 
one  which  is  known  to  be  +1.  This  suggests  that  the  problem  for  a 
nearly  circular  contour  also  has  a favorable  distribution  of  eigen- 
values which  makes  a rapidly  converging  iteration  process  feasible. 

Speaking  in  more  physical  terms,  one  can  say  that  by  mapping 
a profile  onto  a contour  which  is  nearly  circular,  a near  circle, 
one  suppresses  the  interference  between  the  upper  and  lower  sides 
which  is  generated  by  the  form  in  which  the  solution  is  expressed 
(superposition  of  sources  or  vortices).  Therefore,  the  first  step 
of  the  Theodorsen  method,  namely^  transformation  of  the  profile  into 
a near  circle, is  carried  over  in  the  present  work.  The  solution  to 
the  potential  problem  is  then  computed  for  the  near  circle. 

The  benefits  derived  from  the  transformation  into  a near  circle 
are  closely  related.  One  obtains  a contour  which  may  be  considered 
an  analytic  contour.  The  integrand  can  be  represented  by  a very 
smooth  periodic  function.  The  integral  can  then  be  evaluated  very 
accurately  by  one  of  the  simplest  of  all  quadrature  formulas,  namely, 
the  trapezoidal  rule.  The  suppression  of  the  interference  between 
the  lower  and  upper  sides  makes  greater  accuracy  possible  and  reduces 
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the  number  of  quadrature  points  required  to  achieve  a given  accuracy. 
The  resulting  equations  may  be  solved  by  iteration.  Lastly,  the 
Kutta  condition  may  be  applied  directly  to  the  trailing  edge. 

This  raises  the  question  whether  each  element  of  a multielement 
profile  can  be  transformed  into  a near  circle.  The  results  of  [1] 
and  our  own  experience  with  the  Karman  Trefftz-like  mapping,  des- 
cribed in  Appendix  B,  leads  us  to  believe  that  multielement  profiles 
can  be  transformed  into  near  circles  rather  well.  Thus  the  inter- 
ference between  the  upper  and  lower  sides  of  a profile  element  can 
usually  be  suppressed,  but  the  interference  between  portions  of 
elements  close  to  one  another  may  still  cause  difficulties. 

In  many  cases  the  Integral  equation  formulation  is  either  passed 
over  quickly  or  bypassed  by  formulating  an  approximate  expression 
for  the  tangential  velocity  directly  from  physical  considerations. 

The  author  believes  that  it  is  preferable  to  formulate  the  integral 
equation  explicitly  before  approximations  are  Introduced  for  then  all 
approximations  are  made  at  once,  and  it  might  be  possible  to  achieve 
high  accuracy  without  much  additional  effort. 

For  instance,  in  the  present  case  one  must  make  certain  assump- 
tions when  evaluating  numerically  the  Integral  in  the  Integral  equa- 
tion. Frequently,  the  assumption  is  made  that  the  density  is  approx- 
imated by  piecewise  linear  functions.  From  a mathematical  point  of 
view  the  Idea  of  interpolation  by  trigonometric  polynomials  is  attrac- 
tive and  we  use  it  in  the  present  approach.  (As  far  as  the  computa- 
tions are  concerned  this  interpolation  need  not  be  carried  out.) 
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Obviously  this  is  a much  better  representation  of  the  density  func- 
tion than  an  approximation  by  piecewise  linear  functions.  With  this 
approximation  the  integrals  to  be  evaluated  have  very  smooth  inte- 
grands and  the  trapezoidal  rule  gives  much  better  results  than  prac- 
tically all  higher  order  methods.  The  coiuputatlon  to  be  carried  out 
is  extremely  simple,  actually  the  same  as  an  approximation  of  the 
integrands  by  piecewise  linear  functions.  But  the  rationale  given 

I 

shows  that  great  accuracy  can  be  achieved  and  sample  calculations 

[ bear  this  out.  These  refinements  are  perhaps  not  oosslble,  or  in 

} 

I any  case,  not  clearly  recognizeable  if  one  Intrciuces  approximations 

' at  an  earlier  stage. 

-- 

^ If  one  has  a single  profile  then  a conformal  mapping  of  the  pro- 

file on  a circle  is  practical  because  a solution  for  the  potential 
field  around  a circle  is  known.  If  one  proceeds  to  a bielement  air- 
foil then  a mapping  of  the  flow  field  on  the  inside  of  a ring  is 
possible,  in  principle,  and  the  solution  in  this  normalized  form  is 
again  known,  but  only  in  principle.  The  procedure  for  the  bielement 
profile  is  closely  connected  with  the  use  of  elliptic  functions  which 
are  not  always  available  in  a computer.  If  an  airfoil  has  more  than 
two  elements  then  the  idea  of  conformal  mapping  onto  a standard  region 
breaks  down  completely. 

« 

In  contrast  the  Integral  equation  method  can  be  readily  extended 
to  a multielement  section.  It  must  be  admitted  that  for  the  integral 
equation  method  one  must  solve  the  integral  equation  several  times 
to  allow  for  angle  of  attack  and  circulation.  However,  the  kernels 
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of  the  integral  equations  do  not  change  In  this  procedure,  only  the 
right  hand  sides. 

In  the  first  example  treated,  the  author  started  with  the  co- 
ordinates for  profile  points  as  given  in  the  paper  by  Theodorsen 
[2].  These  profile  points  were  transformed  into  a near  circle.  To 
satisfy  the  condition  that  the  normal  velocity  vanish,  one  must  ob- 
tain the  direction  of  the  normal  to  the  near  circle.  For  this  pur- 
pose an  interpolation  by  cubic  splines  of  the  nearly  circular  con- 
tour was  used.  The  graph  of  the  derivatives  of  the  spline  inter- 
polating polynomials  was  ragged.  This  raggedness  is  reflected  in 
the  pressure  coefficient.  One  might  say  that  the  smoothness  of  the 
pressure  coefficient  reflects  the  smoothness  (or  lack  thereof)  of 
the  curve  of  tangential  directions. 

If  the  profile  is  given  merely  by  its  coordinates,  then  one  deals 
with  an  ill-posed  problem,  for  the  steps  which  lie  between  the  data 
and  the  resulting  pressures  Include  a differentiation.  This  means 
that  all  high  accuracy  solutions  of  the  problem  are  bound  to  give  poor 
and  ragged  pressure  coefficients  whenever  the  initial  data  are  poor. 
The  author  rejected  the  idea  of  a posteriori  smoothing  of  the  pressure 
coefficient  in  favor  of  a step  in  which  the  original  profile  data  are 
smoothed.  Such  a step  is  a necessary  preliminary  to  the  determination 
of  the  pressure  coefficient. 

In  this  paper  some  smoothing  has  been  carried  out  in  the  near 
circle  plane  and  after  the  smoothing  the  profile  coordinates  have 
been  recomputed  so  that  one  knows  exactly  which  problem  has  finally 
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been  solved.  Actually,  the  choice  of  the  profile  and  the  smoothing 
procedure  is  not  a mathematical,  but  rather  an  engineering  decision. 
The  graph  of  the  pressure  coefficient.  Figure  1,  Illustrates  what 
can  be  achieved  with  a very  simple  smoothing  process. 

We  have  noted  above  the  benefits  resulting  from  transforming 
profile  elements  into  near  circles.  The  use  of  such  preliminary 
mappings  is  restricted  to  the  two  dimensional  case.  Thus  the  meth- 
ods proposed  here  for  determining  the  pressure  coefficient  on  air- 
foil profile  elements  are  limited  tj  the  two  dimensional  case. 
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SECTION  II 


THE  INTEGRAL  EQUATIONS 

From  a practical  point  of  view  the  primary  objective  is  the 
pressure  distribution  on  an  airfoil  profile.  The  pressure  co- 
efficient can  be  obtained  from  Bernoulli's  equation  when  the  vel- 
ocity 0.1  the  profile  is  known.  In  this  section  we  obtain  expressions 
for  the  theoretical  velocity  at  a point  on  the  profile.  We  examine 
four  possible  ways  for  developing  expressions  for  this  theoretical 
velocity.  It  is  shown  that  there  are  two  essentially  different  ex- 
pressions for  the  velocity. 

The  theoretical  velocity  is  the  solution  of  an  Integral  equation 
or  depends  upon  the  solution  to  an  integral  equation.  In  either 
event  the  properties  of  the  integral  equations  and  their  kernels  are 
well  established  and  can  be  used  to  advantage  for  the  determination 
of  the  numerical  values  of  the  velocity  at  the  profile.  The  two  ex- 
pressions for  the  theoretical  velocity  are  examined  for  possible  ad- 
vantages and  disadvantages. 

The  notation  and  conventions  used  throughout  this  paper  are 
given  in  this  section.  The  pertinent  properties  of  logarithmic  po- 
tential functions  used  in  the  following  developments  are  given  in 
Appendix  A.  Proofs  and  additional  details  on  logarithmic  potential 
functions  may  be  found  in  [3,4,5].  The  reader  who  is  acquainted 
with  the  integral  equation  methods  for  determining  the  flow  field 
around  a body  or  bodies  may  skip  this  section. 


9 


JTECEmiO 


PATE  PLArJK-NoT  FlfJfED 


k- 


Profile  contours  are  usually  quite  slender,  at  least  in  the 
vicinity  of  the  trailing  edge  where  there  is  a sharp  corner.  In 
order  to  obtain  better  conditioned  integral  equations  which  may  be 
solved  by  Iteritlon  the  profile  contour  is  mapped  onto  a nearly  cir- 
cular contour.  The  equation  for  this  nearly  circular  contour  may 
be  considered  analytic.  The  mapping  function  is  conformal  at  all 
points  on  and  exterior  to  the  profile  except  the  trailing  edge. 

Such  a mapping  functflon  is  described  in  Appendix  B.  If  the  airfoil 
profile  consists  of  several  elements,  each  element  is  mapped  into  a 
near  circle  by  repeated  application  of  the  mapping  function.  The 
ability  to  solve  Iteratively  is  not  crucial  for  a single  element 
profile  but  could  be  advantageous  for  multielement  profiles  where  the 
number  of  equations  and  unknowns  may  be  quite  large. 

Let  denote  a simple  closed  and  nearly  circular  contour.  The 
coordinates  of  a point  on  are  given  as  functions  of  the  length  of 
arc  t measured  from  some  convenient  point  on  C^.  (There  should  be 
no  difficulty  in  writing  simply  t rather  than  t^  for  the  length  of 
arc  on  the  i*"*^  contour  C^.)  Hence  let  the  equations  for  be 

q - 

Hf  = n^(t)  for  0 _<  t £ 

where  denotes  the  total  length  of  C^.  It  is  assumed  that  the 
point  ( t ) ,n j (t ) ) traces  in  a counterclockwise  direction  with 
increasing  t.  The  unit  tangent  vector  to  has  components 
(d^^/dt  ,dq^/dt)  and  the  unit  outer  normal  has  components  (drij^/dt. 
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-d^^/dt) . The  functions  5^(t)  and  n^(t)  are  smooth  periodic  func- 
tions of  t.  Lastly,  let  (x^,  y^)  denote  the  coordinates  (or  approx- 
imations to  the  coordinates)  of  the  centroid  of  C^. 

Let  r denote  the  distance  of  the  point  (x,  y)  from  the  point 
(Cj^(t),  n^(t)),  r - /[(x  - + (y  - A real  valued  function 

f(r)  is  then  a function  of  the  four  variables  x,  y,  C,  n*  In  forming 
directional  derivatives  of  f(r)  it  is  helpful  if  one  knows,  without 
depending  upon  context,  which  quantities  are  regarded  as  constants. 
Thus  we  shall  write  for  the  derivative  of  f(r)  in  the  direction  n, 
3f/3nj^  when  (5,  n)  is  fixed  and  3f/3n2  when  (x,  y)  is  fixed. 

The  velocity  field  could  be  determined  if  we  had  either  the  vel- 
ocity potential  function  4i(x,  y)  or  the  stream  function  f(x,  y) . 

We  consider  first  the  determination  of  i^Cx,  y) . The  properties  which 
characterize  i!>(x,  y)  are;  ♦(x,  y)  is  harmonic  outside  C^,  the  normal 

derivative  vanishes  on  and  the  derivatives  34i/3x  and  34>/3y  tend  to 

2 2 

well  defined  limits  as  the  quantity  x + y increases  without  bound. 

It  is  natural  to  represent  ^Cx,  y)  as  a linear  combination  of 
functions  which  have  properties  desired  for  4>(x,  y) . Thus  we  may 
write  ♦(x,  y)  either  as 

y)  = 0(x,  y)  + ?/c.  k(t)[log  r]  dt  (1) 

1 i 
or 

♦ (x,  y)  = I<i(x,  y)  +£  fp  m(t)(3[log  rl/3n„)  dt  (2) 

i ^i  ^ 

where  i goes  from  1 to  N , N is  the  number  of  contours  C..  Here 

c c 1 

i>(x,  y)  denotes  a harmonic  function  which  provides  certain  properties 
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desired  for  it(x,  y)  and  the  Integral  provides  the  capability  for 
satisfying  the  boundary  conditions.  In  Eq . (1)  the  integral  rep- 
resents the  contribution  to  tCx,  y)  of  a continuous  layer  of  sources 
(or  sinks)  of  strength  k(t)  distributed  along  the  contour  C^,  while 
in  Eq.  (2)  the  integral  represents  the  contribution  due  to  a con- 
tinuous layer  of  doublets  of  strength  or  moment  m(t). 

We  consider  first  the  representation  of  ♦(x,  y)  as  given  by 
Eq.  (1).  We  have 


34>(x,  y)/3n  - 84,(x,y)/3n  + J f k(t)(3[log  r]/3n.)  dt  (3) 

i '^i  ^ 

Let  Sj  satisfy  the  condition  0 ^ s^  - ^ j ' point  x * 

y “ -lies  on  the  contour  . We  will  write  <t'(Sj)  for  the  value 

of  ^(x,  y)  at  the  point  (^j(Sj),  nj(s^)).  This  convention  will  be 
used  for  denoting  the  value  of  a function  at  a particular  point  on 
a particular  contour. 

We  obtain  from  Eq . (3)  by  Eq . (A12) 


-34>(s  )/3n  = Tik(s  ) +Zj^  k(t)(3[log  r]/3n  ) dt 
J J i i 


(4) 


since  the  34>^(Sj)/3n  = 0.  The  Eq.  (A)  is  a Fredholm  integral  equation 
of  the  second  kind  with  inhomogeneous  term,  -3^(Sj)/3n.  The  Eq.  (4) 
must  be  solved  for  the  density  k(t)  on  the  contour  C^.  The  Integrands 
are  well  defined  and  continuous  except  for  the  Integral  around  Cj , 
where  the  Integrand  is  Indeterminate  when  t = Sj . The  value  of  this 
integrand  at  t = s^  is  given  in  Appendix  C. 

Once  the  density  k(t)  is  determined  the  velocity  potential  func- 
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r 


t ion  y)  is  known.  The  velocity  at  a point  (^j(Sj), 

the  contour  is  given  by  34>(s^)/aT,  the  value  of  the  derivative 
of  4i(x,  y)  in  the  direction  of  the  tangent  vector  to  at  (^^(Sj), 
Hj (Sj ) ) . We  have 


at 


(s  )/3t  = a.).(s  )/3t  + ZC  k(t)(3[log  rl/ar.)  dt  (5) 
J J i '■ 


Thus  the  evaluation  of  a line  integral  is  required  for  each  poin.  at 
which  the  value  of  the  velocity  is  desired. 

The  integrands  occurring  in  Eq . (5)  are  well  defined  and  contin- 
uous except  for  the  integral  around  the  contour  , where  the  inte- 
grand has  a first  order  singularity  for  t = s^ . In  this  case  one 
must  take  for  the  value  of  the  integral  its  Cauchy  principal  value. 

In  the  numerical  evaluation  of  the  integral  around  the  contour  Cj  we 
use  the  derivative  of  the  density  k(t)  with  respect  to  the  integration 
variable.  The  determination  of  these  derivatives  is  discussed  in 
Appendix  D.  The  numerical  evaluation  of  the  singular  Integral  occur- 
ring in  Eq.  (5)  is  discussed  in  Appendix  E. 

For  our  particular  problem  we  take 


<(i(x,  y)  = Ux  + Vy  + Z Wj  arctan[(y  - y ^)  / (x  - (6) 

where  I runs  through  the  values  1,  2,  ...  , N . The  first  two  terms 

c 

Ux  + Vy  enable  't'(x,  y)  to  satisfy  the  required  conditions  at  a great 
distance  from  the  contours  C^.  The  terms  following  the  summation 
sign  provide  the  circulation.  We  write  the  density  k(t)  as  a sum  also 


k(t)  - -Uk^(t)  + Vk^Ct)  + Z (7) 
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w. 


Wlien  the  value  of  8())(Sj)/3n  as  obtained  from  Eq.  (6)  and  the  value 
of  the  density,  Eq.  (7)  are  substituted  In  Eq.  (A)  one  obtains 
+ 2 integral  equations  for  the  densities  kj^(t),  1^2(1)  and 
^3^(1),  f - 1 N^. 

It  is  clear  that  these  densities  need  to  be  determined  only 
once.  Then  for  any  prescribed  undisturbed  flow  at  infinity  (or  any 
angle  of  attack)  and  prescribed  velocity  at  a point  on  each 
appropriate  values  for  the  coefficients  U,  V and  8,  = 1,  ...  , N^, 

will  give  a potential  function  l'(x,  y),  Eq.  (1)  which  satisfies  the 
prescribed  conditions.  In  particular  the  coefficients  can  be  de- 
termined to  enforce  the  Kutta  condition  at  the  image  of  the  trailing 
edge  on  each  contour  C^.  The  are  determined  by  a system  of  linear 
equations.  The  coefficients  U and  V must  take  into  account  the  dis- 
tortion at  infinity  (if  any)  arising  from  mapping  the  profile  elements 
into  near  circles. 

Now  consider  Eq.  (2), 

'f(x,  y)  = (()(x,  y)  +Z/c  m(t)(3[log  r]/3n2)  dt 

Suppose  i))(x,  y)  harmonic  at  points  (x,  y)  on  and  Inside  the  contour 
1 = 1,  ...  , N^.  It  follows,  since  then  4'(x,  y)  is  harmonic  in- 
side and  341/ 3n  = 0 by  hypothesis  that  4>(x,  y)  is  constant  inside 

Hence  the  limit  of  (x,  y)  as  the  point  (x,  y)  tends  to 

^ constant.  Then  by  Eq . (A8)  we  have 

♦pCSj)  • ♦(Sj)  - irm(Sj)  + X fc  m(t)(3[log  r]/3n2)  dt  (8) 
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By  Eq.  (A6) 


m(t)  » (constant  - 4i^(t))/27t  (9) 

Substitute  for  ra(t)  in  Eq . (8)  its  value  as  given  by  Eq.  (9).  Then 
using  Eq.  (A5)  and  the  fact  that 

J (8[log  r]/3n„)  dt  = 0 if  (x,  y)  lies  outside  C. 

^1 

= Ti  if  (x,  y)  lies  on  C^. 


Eq . (8)  simplifies  to 


2ii(ti(Sj)  = K4>^(Sj)  $^(t)[8B/3t] 


dt 


(10) 


The  Integral  equation,  Eq.  (10),  could  be  solved  for  ♦g(t)  and 
in  this  way  one  could  determine  the  potential  function  $(x,  y) . In- 
stead, differentiate  Eq.  (10)  with  respect  to  . After  an  inte- 
gration by  parts  we  have 


2tii(i's(s  ) = 7i34i  /3s  -XX'  (^^  /8t)(3[log  r]/3n  ) dt 
J ^ ^ i ^ ^ 


(11) 


The  quantity  |34i^/3Sj|  is  the  speed  of  the  flow  past  the  body  bounded 
by  the  contours  C^.  Set  v(t)  » 3<t'^/3t,  then  Eq.  (11)  can  be  lawritten 
as 


2tk|i'  (s  ) = 1TV(S  ) - u v(t)(3[log  r]/3n  ) dt  (12) 

J J i ^i  *■ 

The  Eq.  (12)  is  of  the  same  form  as  Eq.  (All).  From  properties 

of  harmonic  functions  we  conclude  that  <f'(s)  satisfies  the  condition 

expressed  by  Eq.  (A13).  Hence  Eq.  (12)  is  solvable  for  v(t);  however. 
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the  solution  v(t)  is  not  unique.  To  any  solution  of  Eq . (12)  one 
may  add  a solution  of  the  homogeneous  equation.  The  Kutta  condition 
is  used  to  determine  the  solution  of  the  homogeneous  equation  to  be 
added  to  the  solution  of  Eq . (12). 

In  the  method  associated  with  Eq.  (1)  one  provides  for  circu- 
lation by  means  of  the  functions  ui^^arctan  [ (y  - y^)/(x  - 
basic  assumptions  for  the  derivation  of  the  methods  associated  with 
Eq.  (2)  did  not  permit  the  use  of  these  functions.  By  contrast,  the 
circulation  is  provided  for  in  a natural  way  by  the  eigenfunctions 
associated  with  the  eigenvalue  tt.  We  postpone  further  comparison 
until  the  possibilities  for  other  method  are  exhausted. 

The  stream  function  f(x,  y)  is  characterized  by  the  properties: 
y(x,  y)  is  harmonic  at  points  (x,  y)  exterior  to  the  contours  C^, 

y(x,  y)  is  constant  on  and  the  first  order  derivatives  tend  to 

2 2 

finite  limits  as  the  quantity  x + y increases  without  bound.  We 
may  attempt  to  represent  f(x,  y)  in  either  the  form 

'f(x,  y)  = 4i(x,  y)  + Ifr  m(t)(9[log  r]/8n,)  dt  (13) 
1 ^1  ^ 

or 

f(x,  y)  » t|((x,  y)  +Xfr  k(t)[log  r]  dt  (14) 

1 ^1 

The  tangential  velocity  to  the  contours  is  not  readily  computed 
from  f(x,  v)  as  given  by  Eq.  (13).  For  this  reason  we  do  not  con- 
sider the  representation  of  'f(x,  y)  which  contains  the  logarithmic 
potential  of  a double  layer  further. 
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Lastly,  consider  the  stream  function  as  given  hy  Eq.  (14). 
Since  t(x,  y)  is  constant  along  the  contour  C^,  the  derivative  of 
y(x,  y)  at  a point  on  in  the  direction  of  the  tangent  vector  to 
vanishes.  Thus  we  have 


ay(Sj)/3Sj  = 0 = 3il/(Sj)/9Sj  k(t)(3[log  r]/3s^) 


dt 


(15) 


This  is  an  integral  equation  of  the  first  kind  for  k(t).  As  we  have 
already  noted  the  integrand  for  the  line  Integral  around  the 
contour  Cj  has  a first  order  singularity  for  t = s^ . The  Eq.  (15) 
is  not  readily  solved. 

The  condition  on  the  function  values  of  '('(x,  y)  can  be  exchanged 
for  a condition  on  the  normal  derivative  of  y(x,  y)  on  If  i()(x,  y) 

is  harmonic  inside  the  contours  and  continuous  across  the  contours 
then  y(x,  y)  as  defined  by  Eq.  (14)  is  harmonic  inside  the  also. 
Since  y(x,  y)  is  harmonic  inside  the  and  constant  on  the  it 
follows  from  properties  of  harmonic  functions  that  T(x,  y)  is  con- 
stant inside  the  contours  C^.  Hence  3y^(Sj)/3n  = 0 and  we  have  by 
Eq.  (AlO) 


k(t)  - (l/2TT)3y  (t)/3n  (16) 

e 

replacing  the  density  function  by  its  value  as  given  by  Eq.  (16)  we 
have  from  Eq • (A12) 


3'*'e(Sj)/3n  - 34»(Sj)/3n  + (1/2) 3y^(s^ ) /3n  + 

(1/2it)  Z/p  (3y^(t)/3n)(3[log  r]/3nj^)  dt 
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Setting  3'l'^(t)/8n  = -v(t)  and  combining  like  terms  we  obtain 

2nOii((s  )/3n)  = -nv(s  ) v(t)(3[log  r]/3n  ) dt  (17) 

1 1 i 

These  same  equations  for  determining  the  velocity  at  the  contour 
can  be  derived  in  yet  another  way  from  more  elementary  considerations. 
Namely,  one  can  go  from  the  knowledge  of  the  velocity  components  due 
to  a point  source  or  Induced  by  a point  vortex  at  the  origin  to  the 
velocity  components  due  to  a line  distribution  of  sources  or  Induced 
by  a line  distribution  of  vortices.  Usually,  however,  in  such  an 
approach  certain  approximations  are  made  and  accepted  at  an  early 
stage  and,  hence,  the  equations  given  above  are  not  explicitly  ob- 
tained . 

The  Eqs.  (12)  and  (17)  are  essentially  the  same.  Hence  previous 
remarks  concerning  Eq.  (12)  apply  equally  well  to  Eq.  (17).  The  der- 
ivation of  the  equation  for  the  tangential  velocity  starting  with 
Eq.  (2)  and  concluding  with  Eq.  (12)  is  often  referred  to  as  the  meth- 
od of  doublets.  This  last  derivation  starting  with  Eq.  (14)  and 
leading  to  Eq . (17)  is  called  the  method  of  vortices  or  Prager's 
method  [6].  We  will  not  consider  the  method  of  doublets  further. 

The  first  method  which  gave  Eq.  (4)  for  the  density  function  and  Eq. 
(5)  for  the  tangential  velocity  we  will  call  the  method  of  sources. 

In  the  preceding  pages  we  considered  two  representations  for  the 
velocity  potential  and  two  representations  for  the  stream  function. 
These  four  different  starting  points  led  to  two  different  equations 
or  expressions  for  the  velocity  at  the  contour.  Now  we  want  to  con- 


sider  the  determination  of  the  tangential  velocities  further.  Our 
object  here  is  to  obtain  an  estimate  of  work  involved  using  Eq.  (17) 
and  compare  with  a corresponding  estimate  of  work  using  Eqs.  (4)  and 
(5). 

A first  impression  is  that  the  method  of  vortices  is  far  super- 
ior to  the  method  of  sources,  since  the  tangential  velocity  is  ob- 
tained directly  as  the  solution  of  the  integral  equation,  Eq.  (17). 
For  the  method  of  sources,  on  the  other  hand,  the  solution  of  the 
integral  equation,  Eq.  (4)  is  just  the  first  step.  An  integration 
is  required  to  obtain  the  tangential  velocity. 

For  the  method  of  vortices,  let  us  .lote  that  along  with  Eq.  (17) 
we  need  to  consider  the  homogeneous  equation, 

-^V3^(Sj)  V3j^(t)  (3[log  rj/an^)  dt  = 0 (18) 

We  need  independent  solutions  V3^(t),  S,  = 1,  ...  , , of  Eq. 

(18)  in  order  to  be  able  to  enforce  the  Kutta  condition  at  prescribed 
points  on  the  contours  C^.  Also  in  order  to  accomodate  any  given 
undisturbed  flow  at  infinity  write  v(t)  in  Eq.  (17)  as 

v(t)  - Uv3(t)  + Vv2(t)  (19) 

Then  with  i()(x,  y)  = Uy  - Vx  one  obtains  from  Eq.  (17)  one  equation  of 
the  same  form  as  Eq.  (17)  for  v^(t)  and  also  for  V2(t). 

Hence  for  the  method  of  vortices  the  integral  equation,  Eq.  (17) 
has  to  be  solved  for  two  different  inhomogeneous  terms.  The  homo- 
geneous equation,  Eq.  (18)  is  solved  by  iteration.  Note  that  an 
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iteration  is  an  integration.  The  number  of  iterates  depends  upon 
the  tolerance  which  one  sets.  The  total  number  of  integrations  re- 
quired is  then  the  sum  of  the  number  of  iterates  for  the  determination 
of  the  eigenfunctions,  that  is,  the  solutions  of  the  homogeneous 
equation,  F.q . (18). 

For  tne  method  of  sources  the  Integral  equation,  Eq.  (4),  has 
to  be  solved  for  + 2 different  Inhomogeneous  terms.  Next,  N^  + 2 
integrations  have  to  be  performed.  in  the  method  of  vortices  the 
integral  equation  does  not  have  to  be  solved  as  many  times  as  in  the 
method  of  sources.  This  decided  advantage  is  offset  or  diminished 
by  the  number  of  Iterates  required  for  the  determination  of  the 
eigenfunctions. 

We  conclude  that  the  method  of  vortices,  in  general,  requires 
less  computer  resources  than  the  method  of  sources.  Each  method  was 
programmed  so  that  the  system  of  simultaneous  equations  obtained  could 
be  solved  either  directly  or  iteratively.  When  the  number  of  contours 
N(.  is  small  and  high  accuracy  is  desired  the  difference  in  cost  is  neg- 
ligible. As  the  number  of  contours  is  increased  we  expect  the  method  of 
vortices  based  on  Iterative  procedures  to  exhibit  a cost  advantage. 

It  is  clear  that  many  different  computer  programs  for  the  speed  of 
the  flow  past  the  profiles  are  possible  based  on  Eqs . (4)  and  (5)  or  Eqs . 
(17)  and  (18),  respectively.  The  different  programs  are  the  result  of 
different  ways  of  numerically  approximating  and  solving  these  equations. 
Four  programs  are  listed  in  Appendix  G.  These  programs  are  based  on  the 
numerical  approximations  and  solution  methods  discussed  in  Section  III. 
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SECTION  III 


NUMERICAL  APPROXIMATIONS 


In  this  section  we  describe  the  numerical  procedures  we  use  to 
obtain  the  tangential  velocities.  Then  we  mention  briefly  several 
alternative  numerical  procedures.  In  the  previous  section  we  found 
that  whether  one  used  a method  based  on  sources  or  doublets  or  vor- 
tices an  integral  equation  needed  to  be  solved  first.  Note  that  this 
integral  equation,  before  any  modification,  has  the  same  kernel,  re- 
gardless of  method.  Accordingly,  different  numerical  methods  for 
expressing  the  value  of  the  integral  account  for  some  of  the  differ- 
ences in  solution  procedure. 

We  assume  that  the  contours  C^  are  very  smooth.  The  expressions 
for  the  equations  for  the  C^  may  be  considered  analytic.  Then  the 
Integrands  of  the  integrals,  Eqs.  (4),  (17)  and  (18)  are  smooth  con- 
tinuous periodic  functions  of  t for  each  value  of  s^ . These  inte- 
grals, which  are  of  the  form 


/ k(t)(3[log  r]/3nj^)  dt 


are  evaluated  with  high  accuracy  using  the  trapezodial  rule  and 
equally  spaced  Increments  of  the  integration  variable  [7]. 

Let  f(6)  be  a periodic  function  of  6 of  period  2Tt.  Set 
A9  “ 2n/N,  N an  integer  ^ 1 
e(I)  - (I  - 1)A0  for  I = 1,  ...  , N+1 

and 

f(I)  - f(e(D). 
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Then  the  value  of  the  Integral  of  f(6)  over  one  perljd  Is  given 
approximately  by  the  trapezoidal  rule  by  the  expression 
N 

I f(I)A0  (20) 

1=1 

Now  let  f(t)  denote  a function  defined  on  the  simple  closed  con- 
tours C^.  We  found  it  convenient  to  express  each  contour  in 
polar  coordinates  with  pole  at  (or  near)  the  centroid  of  C^.  Then 
with  &0(T)  constant  on  but  not  necessarily  the  sane  constant  on 
different  contours  and  with  appropriate  indexing  one  has  the  approx- 
imation 

NT 

5;]"  f(t)  dt  = I f(i)idt(i)/d0]A0(i)  (21) 


where  NT  denotes  the  total  number  of  stations  on  all  the  contours 
C^.  Thus  the  expression  for  the  numerical  value  of  a line  Integral 
around  one  or  several  contours  is  basically  the  same. 

Expressing  the  values  of  the  integrals  in  Eqs.  (4),  (17)  and 


systems  of  linear  equations.  Using  the  conventions  just  introduced 
the  elements  in  the  coefficient  matrix  A(I,  J)  of  the  system  of  equa- 
tions stemming  from  Eq.  (4)  are,  for  J I,  given  by  the  expression 


(C(l)-^(J))(dn(I)/dt)-(n(I)-n(J))(dg(I)/dt  ^dt(J)^  ^22) 

(C(l)  - C(J))^  + (n(I)  - n(J))^ 

and  the  diagonal  elements  are  given  by 

„ + (d^n(I)/d9^)(dC(I)/d0)  - (d^^(I)/de^)(dn(I)/de)  ^23) 

2{(dC(I)/d0)^  + (dn(I)/de)^l 
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I 


The  right  hand  sides  as  determined  from  Eq  (6)  is 


(dn(I)/de)/(dt(I)/d0),  I = 1,  ... 

, NT 

(24) 

(dC(I)/d0)/(dt(I)/d0),  I = 1,  ... 

, NT 

(2S) 

(C(I)  - iu))(dai)/d0)  + (n(l)  - 

y(S,))(dn(I)/d0) 

(26) 

1(5(1)  - ia))^  + (r,(I)  - v(i)) 

^l(dt(I)/d0) 

for  I = 1 NT,  H = 1,  ...  , N . 

c 

The  off  diagonal  elements  for  the  matrix  A(I,  J)  associated 
with  Eq.  (18)  are  given  by  Eq.  (22).  The  diagonal  elements  are  given 
by  Eq.  (23)  with  the  term  ir  deleted.  The  eigenvectors  associated 
with  the  eigenvalue  w are  determined  by  iteration.  For  the  initial 
guess  or  starting  value  we  take  an  eigenvector  determined  from  an 
eigenfunction  corresponding  to  the  eigenvalue  n of  the  adjoint  ker- 
nel . 

The  set  of  functions  defined  on  the  contours  and  which  have 

the  value  1 on  one  of  the  contours  and  0 on  all  other  contours  is  a 

complete  set  of  eigenfunctions  corresponding  to  the  eigenvalue  ti  of 

the  adjoint  kernel.  Accordingly  the  vector  with  components  equal  to 

1 at  the  stations  on  one  contour  and  zero  otherwise  is  taken  as  the 

initial  value  for  the  iteration  process.  For  our  test  problem  the 

maximum  difference  between  the  components  of  successive  iterates  was 
“8 

less  than  10  after  six  to  eight  iterations. 

Denote  the  components  of  the  eigenvectors  just  determined  by 
CC(I,  J),  I = 1,  ...  , NT  and  J • 1,  ...  , N^.  Denote  the  components 
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I 


of  the  eigenvectors  defined  by  the  known  eigenfunctions  of  the  ad- 
joint kernel  by  CA(I,J),  I + 1,  ...  , NT,  J “ 1,  ...  , N . One 

c 

should  note  that  the  eigenvectors  CA(I,  J)  are  not  eigenvectors  of 
the  transpose  of  the  matrix  A(I,  J).  Instead,  they  are  the  eigen- 
vectors of  the  matrix  In  the  linear  system  approximation,  to  the 
homogeneous  Integral  equation  with  the  adjoint  of  the  kernel  In  Eq. 
(18).  This  difference  between  the  eigenfunctions  of  the  adjoint 
Integral  equation  and  the  eigenfunctions  determined  by  the  transpose 
of  the  matrix  A Is  the  result  of  using  the  integration  variable  6. 

Once  the  Integral  equation  is  reduced  to  a system  of  equations 
one  could  forget  the  origin  of  the  system  of  equations.  We  believe 
that  it  is  preferable,  in  general,  to  use  relations  involving  the 
original  equations.  In  the  present  problem  It  is  advantageous  to 
proceed  as  though  we  were  dealing  directly  with  the  Integral  equa- 
tions rather  then  linear  system  approximations,  since  the  eigen- 
functions of  the  adjoint  kernel  are  known  whereas  the  eigenvectors 
of  the  transpose  of  the  matrix  A(I,J)  are  not  known. 

Once  the  eigenvectors  of  the  matrix  A(I,J)  are  determined  we 
may  regard  the  eigenfunctions  for  the  integral  equation,  Eq.  (18) 
as  known.  Since  now  we  know  the  eigenfunctions  corresponding  to  the 
eigenvalue  n for  the  kernel  of  Eq.  (17)  and  its  adjoint  we  can  modify 
Eq.  (17)  in  the  manner  discussed  in  Appendix  A,  Eq.  (A21).  The  en- 
tries for  the  matrix  Aj_j(I,J)  in  the  linear  system  approximation  of 
the  Integral  equation  with  the  modified  kernel  are,  for  I J,  glv«'n 
by  the  expression 
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(27) 


N 

c 

-A(I,  J)  + TI  X cc(l,  K)  CA(I,  K)(dt(J)/de)A6(J) 

K»1 

and  the  diagonal  elements  are  given  by 

N 

c 

-A(I,  I)  + Tt[l  + z CC(I,  K)  CA(I,  K)(dt(I)/de)Ae(I)]  (28) 

K=1 

where  the  A(I,J)  are  the  elements  of  the  matrix  associated  with  the 
Integral  equation,  Eq . (18).  One  is  assured  that  the  matrix  Aj^(I,J) 
is  nonsingular.  Hence  the  associated  linear  system  can  be  solved 
by  direct  means. 

The  linear  system  associated  with  the  matrix  Aj^(I,J)  has  to  be 
solved  for  two  right  hand  sides.  The  components  for  these  vectors 
are 


2n(d^(I)/de)/(dt(I)/d0),  I - 1 NT  (29) 

and  2w(dn(I)/d0)/(dt(I)/d0),  I - 1 NT  (30) 

The  tangential  derivative  at  a station  on  one  of  the  contours  is 

then  given  as  a linear  combination  of  the  eigenfunctions  CC(I,J) 

and  of  the  solutions  corresponding  to  the  right  hand  sides  given  by 

the  expressions  (29)  and  (30). 

The  system  of  linear  equations  obtained  from  the  Integral  equa- 
tions, either  Eq.  (4)  or  Eq.  (17),  is  well  suited  to  solving  by 
Iterative  means.  In  vector  matrix  notation  the  linear  system  is  of 
the  form 

(XI  + A)  X - b (31) 
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where  X is  the  maxlmuin  eigenvalue  of  A.  Note  that  the  term  ttI  is 
incorporated  in  the  matrix  defined  by  Eqs.  (22)  and  (23).  Divide 
Eq.  (31)  by  X and  consider  first  the  case  with  the  plus  sign  in 
front  of  the  matrix  A.  Thus  Eq.  (31)  becomes  (I  + A/X)  x « b/X  and 
a linear  system  of  this  form  represents  the  integral  equation,  Eq. 
(4).  Formally  one  has 

X = b/X  - (A/X)(b/X)  + (A/X)^(b/X)  - + ...  (32) 

If  the  vector  b has  no  component  in  space  spanned  by  the  eigen- 
vectors corresponding  to  the  eigenvalue  X the  right-hand  side  of  Eq. 
(32)  will  converge  to  a solution  of  Eq.  (31).  (We  have  assumed  that 
-X  is  not  an  eigenvalue  of  the  matrix  A.) 

To  Insure  that  the  right-hand  side  has  no  component  in  the  space 
spanned  by  the  eigenvectors  corresponding  to  X write  x as 

x ” X + a(b/X)  (33) 

Then  we  can  exchange  the  system  of  equations  for  the  unknown  x for 
the  following  system  of  equations  for  the  unknown  x 

(1  + A/X)i  ■=  b/X  - a(b/X)  - a(A/X)(b/X)  (34) 

Observe  that  if  a « 1/2  the  right-hand  side  of  Eq.  (34)  has  no  com- 
ponent in  the  space  spanned  by  the  eigenvectors  of  A corresponding 
to  the  eigenvalue  X.  Thus,  replacing  (b/X)  in  Eq.  (32)  by  the  right- 
hand  side  of  Eq.  (34)  gives  x and  then  Eq.  (33)  gives  x. 

The  case  of  the  minus  sign  in  Eq.  (31)  is  similar,  up  to  a 
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point,  to  that  of  the  plus  sign.  The  equation  (I  - A/X)  x = b/A 
represents  the  integral  equation  Eq.  (17).  As  for  the  plus  sign, 
this  equation  has  solutions  only  if  the  right-hand  side  is  orthogonal 
to  every  eigenfunction  corresponding  to  X of  the  adjoint  kernel. 
Hence,  if  the  right-hand  side  b/X  has  no  component  in  the  space 
spanned  by  the  eigenfunctions  of  A corresponding  to  X,  the  expansion 

X = (b/X)  + (A/X)(b/X)  + (A/X)^(b/X)  + ...  (35) 

should  converge. 

If  the  right-hand  side  b/X  has  a component,  in  the  space  of  the 
eigenvectors  of  the  matrix  A,  sufficiently  large  to  affect  the  con- 
vergence to  X in  Eq.  (35)  the  procedure  given  by  Eqs.  (33)  and  (3A) 
cannot  be  used  to  remove  this  component.  Considerable  additional 

computations  are  required.  First  one  would  have  to  obtain  the 

T T 

eigenvectors  of  A . Next,  these  eigenvectors  of  A would  be  used 

to  obtain  a system  of  linear  equations  for  the  component  of  b/X 

in  the  space  spanned  by  the  eigenvectors  of  A.  Then  this  system 

of  equations  would  have  to  be  solved  for  this  component. 

The  tangential  velocities  can  be  obtained  from  either  the  vel- 
ocity potential  function  *(x,y)  or  the  stream  function  4'(x,y).  Using 
a continuous  distribution  of  sources  in  the  representation  for  the 
potential  function  led  to  an  integral  equation  which  had  to  be  solved 
in  order  to  determine  the  tangential  velocities.  Using  a continuous 
distribution  of  doublets  to  represent  the  potential  function  or  a 
continuous  distribution  of  vortices  to  represent  the  stream  function 


led  to  an  essentially  different  Integral  equation.  Thus  there  are 
two  basic  approaches  in  the  determination  of  the  tangential  vel- 
ocities. , 

In  the  above  pages  of  this  section  we  described  two  ways  for 
determining  a numerical  solution  of  each  Integral  equation.  In  all 
cases  the  value  of  the  Integral  was  expressed  by  the  trapezoidal 
rule  with  equal  increments  of  the  integration  variable.  Every  dif- 
ferent method  for  expressing  the  value  of  the  Integral  in  either  in- 
tegral equations  may  be  considered  a different  procedure  for  deter- 
mining the  tangential  velocities. 

Instead  of  using  continuous  distributions  of  sources,  doublets 
or  vortices  one  might,  for  example,  approximate  the  potential  func- 
tion with  finitely  many  sources,  with  strengths  to  be  determined, 
distributed  along  the  contours  C^.  In  order  to  determine  the  unknown 
strengths  the  expression  for  the  normal  velocity  at  points  on  the 
midway  between  two  sources  is  equated  to  zero.  In  this  way  one 
obtains  a system  of  equations  for  the  source  strength.  One  could 
proceed  in  an  analogous  way  using  either  vortices  or  doublets,  [8, 
pp.  193-6].  Using  a discrete  rather  than  a continuous  distribution 
gives  the  appearance  of  avoiding  the  Integral  equation  and  the  prob- 
lem of  evaluating  the  Integral  associated  therewith.  However,  one 
pays  a price  for  this  computational  simplicity.  The  quantities  of 
interest  are  the  velocities  at  points  along  the  contours.  The  lumping 
of  the  singularities  into  discrete  points  on  the  boundaries  ad- 
versely affects  the  accuracy  of  this  method  in  the  vicinity  of  the 
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C^.  It  is  clear  that  at  points  where  the  sources  (vortices  and 
doublets)  are  located  the  velocities  are  infinite. 

Another  method  for  evaluating  numerically  the  line  Integrals 
occurring  in  Eqs.  (4)  and  (17)  is  the  "panel"  method.  In  this  method 
the  contours  are  replaced  by  polygonal  contours  and  the  density 
k(t)  (method  of  sources)  or  the  tangential  velocity  v(t)  (method  of 
vortices)  is  taken  as  constant  on  each  side  of  the  polygonal  con- 
tours. (Different  constant  on  different's'ides,  of  course.)  Instead 
of  constant  functions  on  the  sides  of  the  polygonal  contours  one  can 
use  linear  functions.  Thus  for  the  side  of  length  I between  the 
Ith  and  the  (I  + 1)  ^ stations  take 

k(t)  = k(I)  + [k(I  + 1)  - k(I)Ht/i) 

For  0 £ t ^ i 

The  resulting  line  integrals  along  the  sides  of  the  polygonal  con- 
tours are  well  defined  and  finite  for  all  points  except  the  end 
points  of  the  sides  of  the  polygonal  contours. 

The  details  for  evaluation  of  such  an  integral  with  uniform 
density  are  given  in  [8,  pp.  191-3]  for  the  special  case  of  a seg- 
ment of  the  x-axls.  Using  a rotation  of  the  xy-plane  one  can  re- 
duce the  evaluation  of  the  integral  over  an  arbitrarily  oriented 
side  to  this  special  case. 

Details  (which  Include  the  case  of  a constant  density)  for  the 
case  of  a density  linear  on  a side  are  given  in  [9,  pp.  21-32]. 

One  can  give  expressions  for  the  velocity  components  directly  in 


terms  of  sums  of  lategrals  of  the  form  being  considered  here.  In 
this  way  one  appears  to  have  avoided  the  integral  equations  of  Sec- 
tion II.  In  actuality  one  has  gone  directly  to  an  approximate  sol- 
ution of  one  of  the  Integral  equations.  However,  in  the  process  of 
going  directly  to  an  approximate  solution  it  may  not  be  clear  what 
approximations  were  made.  This  method  of  evaluating  the  integrals 
in  Eqs.  (4)  and  (17)  or  for  deriving  an  approximate  solution  is 
attractive  because  the  values  of  the  Integrals  are  given  in  terms 
of  readily  evaluated  functions.  One  avoids  the  Indeterminate  case 
occurring  in  the  use  of  the  trapezoidal  rule.  The  determination 
of  the  tangential  velocities  by  the  panel  method  Involves  the  eval- 
uation of  essentially  the  same  Integrals  as  in  the  determination 
of  the  densities.  These  attractive  features  are  offset  by  the  fact 
that  this  method  of  evaluating  the  integrals  requires  a rather  large 
number  of  stations  on  the  contours  to  achieve  a given  accuracy. 
Application  of  the  Kutta  condition  requires  a modified  treatment. 
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SECTION  IV 


DISCUSSIONS  AND  CONCLUSIONS 

The  procedures  described  in  Section  III  were  applied  to  the 
data  lor  a single  element  profile  given  by  Theodorsen  [2j.  Figure 
i IS  the  gr.iph  of  the  pressure  coefficient  P using  40  points  around 
the  near  circle  contour.  (All  the  figures  in  this  report  were  ma- 
chine drawn  and  not  retouched.)  The  appearance  of  Figure  1 could 
be  improved  by  computing  more  values  of  the  pressure  coefficient 
in  the  vicinity  of  the  leading  cage. 

The  present  results  and  those  of  Theodorsen  are  fairly  close. 
This  shows  mainly,  that  the  program  is  trustworthy.  As  no  exact 
results  are  available,  one  cannot  say  that  one  method  gives  greater 
accuracy  than  the  other.  Actually  both  methods  can  be  perfected  to 
give  a high  degree  of  accuracy.  The  present  method  is  advocated 
because  it  extends  readily  to  multielement  profiles  rather  than  be- 
cause of  its  inherent  accuracy. 

The  distribution  of  pressures  shown  by  Figure  1 has  been  deter- 
mined for  the  profile.  Figure  2.  This  profile  is  obtained  from  the 
original  profile  coordinates  Figure  3 by  a smoothing  process.  One 
observes  that  the  "smoothed"  profile.  Figure  2,  is  thicker  than  the 
original  profile.  Figure  3.  It  is,  of  course,  clear  to  the  reader 
that  this  thickening  is  a result  of  an  imperfect  smoothing  process 
and  not  of  the  transformation  to  the  near  circle.  A smoothing  pro- 
cess is  needed  because  pressure  results  will  reflect  irregularities 
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in  the  profile  slope.  Such  irregularities  occur  rather  easily  if 
the  profile  is  described  by  its  coordinates. 

In  Figure  4 we  give  graphs  of  the  quantities  R-1  versus  G and 
dR/dG  versus  0,  where  R is  the  near  circle  radius  before  smoothing. 

One  observes  that  the  graph  of  dR/dG  is  rather  wavy.  Theory  asserts 
and  computational  experience  shows  that  a high  accuracy  method  will 
produce  similar  Irregularities  in  the  pressure  coefficient.  Fre- 
quently the  irregularities  occurring  in  the  pressure  coefficient 
are  removed  rather  arbitrarily  either  on  the  drawing  board  or  by 
truncating  the  expression  for  the  pressure  coefficient.  Thus  cruder 
methods  may  appear  better  because  they  are  less  likely  to  pick  up 
irregularities  in  the  profile  data. 

In  Figure  5 we  give  graphs  of  dR/dG  and  smoothed  values  of 
dR/dG.  From  the  smoothed  values  of  dR/dG  we  redetermined  R,  and 
the  profile  shown  in  Figure  2.  From  this  revised  data  we  determined 
the  pressure  coefficient  shown  in  Figure  1.  We  feel  that  the  appear- 
ance of  Figure  1 implies  that  this  is  the  proper  order  of  proceeding. 
We  can  assert  that  the  data  for  the  profile  of  Figure  2 produced  the 
values  for  the  pressure  coefficient  exhibited  in  Figure  1.  Admit- 
tedly, from  a practical  point  of  view,  the  differences  between  smooth- 
ing the  original  data  and  the  arbitrary  smoothing  of  the  final  re- 
sults in  most  cases  will  not  bo  significant. 

Our  objective  in  applying  a smoothing  process  was  a pressure 
coefficient  graph.  Figure  1,  with  no  undue  irregularities  or  wavl- 
nes'. . We  wanted  to  illustrate  what  a preliminary  smoothing  could 


Resorting  to  a smoothing  process  is  an  aamission  of 
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accoiT.p  i isn . 

some  uisatisfaction  with  the  original  data  or  of  the  fitting  pro- 
cess or  both.  Thus  one  would  expect  the  resulting  profile  after 
smoothing  the  derivatives  to  differ  some  from  the  original.  The 
extent  of  the  difference  wnich  is  acceptable  between  the  original 
profile  and  the  profile  after  smoothing  the  derivatives  is  an  en- 
gineering decision.  It  seems  clear  that,  in  general,  one  should  be 
able  to  determine  a smoothing  for  which  the  resulting  profile  is 
acceptable . 

Our  next  objective  is  a test  of  the  accuracy  of  the  methods 
we  proposed  for  solving  the  integral  equations  and  computing  the 
tangential  velocities.  From  the  above  discussion  it  is  clear  that 
a practical  problem  is  not  well  suited  to  this  task.  For  this  pur- 
pose we  considered  the  complex  potential  function 

F(z)  = z + 0.2/(l-z)  + 0.2/(l+z)  (36) 

This  function  describes  a flow  around  two  nearly  circular  obstacles. 
The  flow  at  infinity  is  parallel  to  the  x-axis  and  of  unit  speed. 

A problem  of  this  kind  would  be  encountered  in  the  determination  of 
the  distribution  of  the  pressure  coefficient  on  biplane  profiles. 

The  "profile"  data  was  determined  from  the  equation 

'f(x,y)  = Im  [F(z)]  = 0 (37) 

That  is,  we  determined  80  points,  40  around  each  contour,  which  sat- 
isfied Eq.  (37).  Figure  6 is  the  graph  of  f(x,y)  = 0.  The  "true 
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value?;"  of  the  tangential  velocities  at  these  80  points  can  be  de- 
termined from  either  Re[F(z)]  or  Im[F(z)],  Eq . (36).  The  maximum 
difference  between  the  tangential  velocities  determined  by  the  nu- 
merical procedures  described  in  Section  III  and  the  true  values  was 
less  than  0.5  E-08. 

r'he  maximum  difference  between  the  true  values  and  the  computed 
values  using  only  20  points,  instead  of  40,  around  each  contour  was 
less  than  0.5  E-08  also.  The  maximum  differences  using  20  points 
was  larger  than  the  corresponding  maximum  difference  using  40  points, 
in  general.  Normally,  one  would  expect  that  doubling  the  number  of 
points  would  produce  a greater  decrease  in  the  maximum  difference. 

We  believe  that  limitations  in  the  accuracy  of  the  initial  data 
accounts  for  the  present  results. 

We  have,  given  two  basic  ways  for  computing  tangential  velocities. 
In  each  of  these  procedures  we  had  to  solve  a system  of  linear  equa- 
tions fseveral  times,  in  fact).  Utilizing  some  special  properties 
of  this  system  of  equations  we  are  able  to  solve  this  system  of 
equations  Iteratively  as  well  as  directly.  Thus  we  may  say  that  we 
have  four  ways  for  computing  tangential  velocities.  All  four  ways 
are  equally  accurate.  The  method  based  on  a vortex  distribution 
(Prager's  method)  required  slightly  less  central  processing  time 
than  did  the  method  of  sources.  With  40  points  per  contour  iter- 
ative methods  used  less  central  processing  time  than  direct  methods. 
For  20  points  per  contour  the  direct  methods  used  less  central  pro- 
of ■ s ,s  i r t imr . 


Usin^’,  a source  distribution  and  solving  the  system  of  linear 
equations  directly  took  1.129  sec.  of  central  processing  time.  With 
a vortex  distribution  and  solving  iteratively  1.010  sec.  of  central 
processing  time  was  used.  These  times  are  for  the  case  of  40  points 
per  contour.  This  central  processing  time  difference  is  not  signif- 
icant. We  expect  that  with  a greater  number  of  profile  elements  the 
method  using  a vortex  distribution  and  solving  the  linear  systems 
iteratively  will  improve  its  advantage  over  the  other  three  methods. 
We  note  also  that  solving  a system  of  equations  iteratively  is  more 
flexible  than  solving  the  same  system  directly^  in  the  sense  that 
the  solving  process  can  be  shortened  by  reducing  the  accuracy  re- 
quirement . 

Thus  we  conclude  that  all  four  methods  proposed  for  the  task 
of  providing  the  basic  Integral  equation,  solving  this  equation  and 
determining  the  tangential  velocities  in  a transformed  plane  are 
good  and  very  accurate.  The  method  based  on  a vortex  distribution 
and  which  solves  the  linear  system  of  equations  iteratively  has  a 
slight  cost  advantage. 

As  we  have  mentioned  the  complete  procedure  requires  some  addi- 
tional steps:  First,  the  mapping  of  the  profile  data  from  the  phys- 

ical plane  to  the  plane  of  near  circlesj  secondly,  the  fitting  of 
the  near  circle  data  and  smoothing.  Appendix  F.  As  a result  of  these 
two  steps  one  removes  the  difficulties  associated  with  the  sharp 
trailing  edge  and  the  leading  edge.  One  has  to  deal  with  very  smooth 
functions.  This  simplifies  solving  the  associated  integral  equation 


and  enables  one  to  obtain  a very  accurate  solution  from  relatively 
few  quadrature  points.  The  Kutta  condition  can  be  invoked  at  Che 
trailing  edge  in  a straight  forward  manner. 

The  disadvantages  of  these  two  steps  are  clear.  In  general, 
they  are  not  fully  automated.  In  particular,  the  smoothing  may  re- 
quire some  drawing  and  reading  of  points  with  a digitizer.  On  the 
other  hand,  for  example,  the  panel  methods  which  bypass  these  first 
two  steps  require  considerably  more  data  points,  especially  along 
the  leading  edge.  Some  times  other  special  precautions  have  to  be 
taken  to  avoid  ill  conditioned  or  singular  coefficient  matrices  in 
the  associated  systems  of  linear  equations. 


Q DENOTES  DRIEINRL  DRTR. 

— INVERSE  KRRMRN  TREFFTZ  MRPPING 
OF  THE  FIT  OF  THE  NERR  CIRCLE 
DRTR  AFTER  SMOOTHING. 


Figure  2. 

The  "Smoothed"  Profile 


38 


0 DENOTES  DRIB  I NHL  DHTR. 

INVERSE  KRRMHN  TREFFT2  MAPPING 

OF  THE  CLHSSICHL  SPLINE  FIT  OF 
THE  NEHR  CIRCLE  DHTR. 


Figure  3. 

The  Clark  Y Profile 
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Figure  4. 

Cubic  Spline  Fit  of  the  Near  Circle  Data  and  Its  Derivative 
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Figure  5. 

The  Derivative  of  the  Cubic  Spline  Fit  of  the 
Near  Circle  Data  and  Its  Smoothed  Replacement 
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Figure  6. 

Graph  of  the  Imaginary  Part  of 
z + 0.2/(l  + z)  + 0.2/(l  - z)  = 0 
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Figure  8. 

The  Image  of  the  Clark  Y Profile  Data 

I 

Under  the  Karman  Trefftz-llke  Mapping 
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APPENDIX  A 


LOGARITHMIC  POTENTIALS 

In  this  appendix  we  give  some  properties  of  harmonic  functions 
expressed  as  the  logarithmic  potential  of  a single  or  double  layer. 
These  properties  are  the  basis  for  the  derivatives  presented  In  Sec- 
tion II.  For  additional  details  and  proofs  see  [3,  4,  5].  We  dis- 
cuss also  the  solvability  of  some  integral  equations  in  this  back- 
ground material  and  a modification  of  these  equations  which  removes 
the  solvability  restriction.  Our  notation,  conventions  and  basic 
assumptions  are  the  same  as  in  Section  II. 

The  function 

w(x,y)  = J m(t)0[log  r]/0n„)  dt  (Al) 

C 

defined  by  the  line  Integral  is  called  the  logarithmic  potential  of 
a double  layer  of  density  (or  moment)  m(t)  distributed  over  the  con- 
tour C.  The  density  function  m(t)  is  supposed  continuous.  The  in- 
tegral in  Eq.  (Al)  defines  a function  w(x,y)  which  is  continuous, 
has  continuous  derivatives  of  all  orders  and  satisfies  Laplace's 
equation 

Aw  = 8^w/3x^  + 9^w/3y^  = 0 
at  all  points  (x,y)  not  on  C. 

The  Integrand  in  Eq.  (Al)  may  be  written  in  other  ways,  such 
as 
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m(t) (l/r)3r/Dn2  (A2) 
m(t)((y  - n(t))(d^/dt)  - (x  - ?(t)) (dn/dt))/r^  (A3) 
m(t) (cos(r ,n))/r  (A4) 
m(t)d6/dt  (A5) 


In  Eq.  (A4)  cos(r,n)  denotes  the  cosine  of  the  angle  between  the 
normal  and  the  radius  vector  to  the  point  (C.h)  from  the  point  (x,y). 
In  Eq.  (A5)  f?  is  the  angle  6 - tt  where  6 = arctan(y  - n)/(x  - O- 
The  function  w(x,y)  does  not  vary  in  a continuous  fashion  as 
the  point  (x,y)  passes  through  the  contour  C.  Let  s denote  an  ar- 
bitrary but  fixed  value  satisfying  0 < s < L.  Let  w. (s)  and  w (s) 

— — i e 

denote  the  limits  of  w(x,y)  as  the  point  with  coordinates  (x,y) 
tends  to  the  point  (^(s),  n(s))  on  C through  points  lying  only  in- 
side C and  points  lying  only  outside  C,  respectively.  Then  w^(s), 
Wg(s)  and  m(s),  the  value  of  m(t)  at  (C(s),  n(s))  on  C,  satisfy 
the  following  conditions; 


w^(s)  - w^(s)  • 2itm(s) 

(A6) 

w (s)  = nm(s)  +J  m(t)(3[log  r]/8n„)  dt 

(A7) 

^ C ^ 

W (s)  = -irm(s)  +J  m(t)0[log  r]/3n„)  dt 
C ^ 

(A8) 

at  (?(s),  ri(s)).  Equations  (A6) , (A7)  and  (A8)  exhibit  the  connec- 
tion between  the  density  function  m(t)  and  the  boundary  values  of 
the  function  w(x,y). 

In  the  process  of  deriving  the  relations  (A6),(A7)  and  (A8)  it 
is  seen  that  it  is  an  eigenvalue  belonging  to  the  kernel  3[log  r]/8n2 
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and  that  m(t)  equal  to  a constant  is  an  eigenfunction  corresponding 
to  the  eigenvalue  -n . The  Eqs.  (A7)  and  (A8)  can  be  rewltten  so 
that  the  kernel  Is  (l/w)3[log  r]/3n2.  Then  by  an  argument  similar 
to  that  of  [5,  pp.  309-310]  one  can  show  that  all  the  eigenvalues 
belonging  to  the  "normalized"  kernel  (l/iT)3[log  r]/3n2  are  real  and 
lie  between  -1  and  1.  It  Is  shown  also  in  [4,  p.  104]  that  -1  is 
not  an  eigenvalue  of  the  normalized  kernel. 

Given  the  boundary  values,  either  w^(s)  or  w^(s),  Eqs.  (A7)  and 
(A8)  are  Fredholm  integral  equations  of  the  second  kind  for  the  den- 
sity function  m(t).  From  the  Fredholm  theory  it  is  known  that  Eq . 
(A7)  is  uniquely  solvable  for  any  prescribed  Inhomogeneous  term 
Wj^(s).  The  Eq.  (A8) , on  the  other  hand,  is  not  solvable  unless 
w^(s)  satisfies  an  orthogonality  condition.  Even  when  solvable  the 
solution  is  not  unique,  but  any  two  solutions  differ  at  most  by  a 
constant . 

The  function 

v(x,y)  =/  k(t)  log  r dt  (A9) 

C 

is  called  the  logarithmic  potential  of  a single  layer  of  density 
k(t)  distributed  over  C.  This  function  defined  by  the  line  integral, 
Eq.  (A9),  is  continuous,  has  continuous  derivatives  of  all  orders 
and  satisfies  Laplace's  equation  at  all  points  (x,y)  not  on  C.  The 
function  v(x,y)  is  continuous  even  as  the  point  (x,y)  passes  through 
C,  but  the  normal  derivative  3v/3n^  is  discontinuous  as  the  point 
(x,y)  passes  through  C. 
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Let  s have  a fixed  value  satisfying  the  condition  0 < s < L. 

Let  the  point  (x,y)  be  confined  to  the  normal  line  to  C at  the  point 

(C(s),  n(s)).  Let  3v  (s)/3n  and  3v  (s)/3n,  denote  the  limits  of 

11  el 

3v/3n^  as  the  point  (x,y)  tends  to  (^(s),  n(s))  through  points  (x,y) 
inside  C and  outside  C respectively.  The  values  3Vj^(s) /Sn^^ , 3v^(s)/ 
3uj^  and  k(s)  satisfy  the  following  conditions: 

3v^(s)/3nj^  - 3v^(s)/3n^  = 2nk(s)  (AlO) 

3v.(s)/3n.  = -Tik(s)  +/  k(t)(3[log  r]/3n  )dt  (All) 

c '■ 

3v  (s)/3n  = nk(s)  + J k(t)(3[log  r]/3n  ) dt  (A12) 

C ^ 

The  kernel  3[log  r]/3nj^  is  the  adjoint  of  3[log  r]/3n2.  Hence 
it  has  exactly  the  same  eigenvalues  as  3 [log  r]/3n2.  Given  normal 
derivatives,  either  3v^(s)/3nj^  or  3v^(s)/3nj^,  the  Eqs.  (All)  and 
(A12)  are  Fredholm  integral  equations  of  the  second  kind  for  the 
density  function  k(t) . Equation  (A12)  is  solvable  for  arbitrary 
inhomogeneous  terms  3v^(s)/3n^^.  Equation  (All)  is  solvable  provided 

/ (3v  (t)/3n  ) dt  = 0 (A13) 

C 

and  when  solvable,  the  solution  is  not  unique. 

Observe  that  the  first  order  derivatives  of  w(x,y)  and  v(x,y) 

2 2 

tend  to  zero  as  the  quantity  x + y increases  without  bound. 

Next  we  want  to  make  some  observations  and  remarks  concerning 
the  solutions  and  the  determination  of  solutions  of  the  integral 
equations,  Eqs.  (A8)  and  (All).  When  the  solvability  condition  is 
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satisfied  each  of  these  equations  is  analogous  to  a consistent  but 
dependent  system  of  linear  equations.  As  is  well  known,  such  a sys- 
tem of  linear  equations  does  not  have  an  unique  solution  but  a linear 
manifold  of  solutions.  Occasionally,  one  has  available  information 
which  can  be  used  to  change  the  consistent  dependent  system  of  equa- 
tions into  a consistent  system  of  equations  which  determines  a par- 
ticular solution  in  the  linear  manifold  of  solutions  of  the  depend- 
ent system.  Then  the  linear  manifold  of  solutions  of  the  dependent 
system  is  obtained  by  adjoining  to  the  particular  solution  the  null 
space  of  the  dependent  system. 

It  is  deslreable  to  change  a dependent  system  Into  an  inde- 
pendent system  because  linear  equation  solvers  are  designed  pri- 
marily for  independent  systems.  Frequently  this  change  is  accom- 
plished by  deleting  an  appropriate  number  of  variables  and  equations. 
The  question  then  arises  as  to  which  variables  and  equations  to  de- 
lete as  some  sets  may  not  span  as  well  as  others. 

In  numerical  work  a consistent  dependent  system  may  manifest 
Itself  as  an  independent  system  which  may  be  ill  conditioned.  Even 
though  one  may  get  by  treating  such  a system  directly  [101  we  be- 
lieve a procedure  which  is  correct  in  principle  is  preferable. 

Given  a dependent  but  consistent  system  of  linear  equations  the 
process  of  deleting  variables  and  equations  is  correct  in  principle. 
However,  in  the  present  case  in  which  the  linear  system  is  obtained 
from  an  Integral  equation  we  believe  an  alternative  procedure  to  be 
described  next  is  preferable. 
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We  have  noted  that  the  kernels  of  Eqs.  (AS)  and  (All)  are  ad- 
joint to  one  another,  it  is  a simple  eigenvalue  belonging  to  these 
kernels  and  for  any  constant  c/0,  f(t)  = c is  an  eigenfunction  or 
solution  of  the  homogeneous  equation  Eq.  (AS).  Usually  we  take 
c * 1.  Let  g(t)  denote  an  eigenfunction  or  solution  of  the  homo- 
geneous equation  Eq.  (All).  The  function  g(t)  is  not  explicitly 
known. 

From  the  Fredholm  theory  we  know  that  for  any  continuous  func- 
tion k(t)  the  image  v(s)  where 

v(s)  = Tik(s)  -/  k(t)(3[log  r]/3n^)  dt  (A14) 

C 

satisfies  the  orthogonality  condition 


f v(s)  f(s)  ds  =J  v(s)  ds  = 0 
C C 

and 

f g(s)  f(s)  ds  =f  g(s)  ds  # 0 
•^C  C 

Consider  now  the  Fredholm  integral  equation 


(A15) 

(A16) 


u(s)  = TTk(s)  -/  k(t)(3[log  r]/3n  ) dt  - itj  k(t)  dt  (A17) 
C ^ C 

= v(s)  - T'J  k(t)  dt 
C 

According  to  the  Fredholm  alternative  either  Eq.  (A17)  has  a unique 
solution  for  any  inhomogeneous  term  u(s)  or  the  homogeneous  equation 
has  non  trivial  solutions.  It  is  readily  seen  that  assuming  the 
homogeneous  equation  has  non  trivial  solutions  leads  to  absurdities. 
Hence  it  follows  Eq.  (A17)  has  a unique  solution  for  every  inhomo- 
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geneous  term  u(s). 

Next  suppose  the  inhomogeneous  term  u(s)  in  Eq , (A17)  satisfies 
the  orthogonality  condition 

J u(t)  dt  = 0 

C 

Then  the  solution  k(s)  of  Eq.  (A17)  satisfies  this  orthogonality 
condition  also.  (Assuming  otherwise  gives  rise  to  an  absurdity  as 
above.)  Hence  k(t)  satisfies  the  integral  equation,  Eq.  (All)  and 
the  manifold  of  solutions  of  Eq.  (All)  is  given  by 

fe(s)  = k(s)  + const.  g(s)  (A18) 

Since  Eq.  (A8)  does  not  arise  in  the  investigations  carried 
out  in  the  remainder  of  this  paper  we  could  omit  further  consider- 
ation here.  However,  for  completeness  and  to  increase  understanding 
of  the  principles  involved  we  shall  discuss  Eq.  (A8)  and  its  modi- 
fication briefly.  The  modification  appears  to  be  the  same  as  for 
Eq.  (All);  that  is,  we  include  in  the  right-hand  side  of  Eq.  (A8) 
a term  of  the  form 

- tj  m(t)  dt 
C 

The  modification  in  both  instances  consists  of  adding  to  the  or- 
iginal kernel  a degenerate  kernel  of  rank  1.  This  fact  is  hidden 
or  may  be  overlooked  because  the  degenerate  kernel  is  the  product 
of  two  constant  functions.  Vfhile  the  degenerate  kernel  is  the  same 
in  both  cases,  its  relation  to  the  original  kernels  differs.  In  the 
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one  case  the  degenerate  kernel  consists  of  the  self  product  of  the 
eigenfunction  f(t)  = 1 of  the  adjoint  kernel  while  in  the  other  case 
it  is  the  self  product  of  the  eigenfunction  f(t)  = 1 of  the  kernel 
itself.  The  factors  in  the  degenerate  kernel  are  determined  by  the 
information  available. 

The  modification  term  is  of  the  form 


J f,(s)  f (t)  k(t)  dt  (A20) 

C 

Ideally,  the  function  fj^(t)  is  an  eigenfunction  of  the  kernel  being 
modified  and  f2(t)  is  an  eigenfunction  of  the  adjoint  kernel.  (It 
is  understood  that  f^^  and  f^  correspond  to  the  same  eigenvalue.) 

In  general  the  desired  eigenfunctions  are  not  known.  This  modifi- 
cation technique  can  still  be  used  if  one  knows  functions  made  up 
from  one  or  the  othev  of  the  desired  eigenfunctions. 

If  the  eigenvalue  is  of  multiplicity  p then  the  modification 
consists  of  a sum  of  p terms 


I./.  ff.Cs)  f2,(t)  k(t)  dt 


j = l C 


The  functions  fj^j(t),  j = l,  ...  , p must  be  linearly  independent  and 
so  also  the  functions  f2j(t). 


APPENDIX  B 


THE  KARMAN  TREFFTZ-LIKE  TRANSFORMATION 

In  this  section  we  display  the  Karman  Trefftz-llke  transfor- 
mation. We  noted  that  the  first  step  in  our  procedure  is  the  con- 
formal transformation  from  the  physical  plane  to  the  near  circle 
plane.  The  purpose  of  this  transformation  is  to  remove  the  sharp 
trailing  edge  and  to  reduce  the  interference  between  the  sources 
over  those  portions  of  the  profile  where  the  upper  and  lower  sur- 
faces are  close  together.  Once  this  transformation  into  a near 
circle  has  been  achieved  no  further  transforming  is  carried  out. 

By  contrast,  in  the  Theodorsen  method  [2]  a transformation  into  a 
near  circle  is  a first  step  in  the  transformation  into  a circle. 

The  Karman  Trefftz  transformation  is  a familiar  tool  of  hydro- 
dynamics. In  this  section  it  is  broken  down  into  a sequence  of 
elementary  transformations.  In  this  form  one  can  recognize  rather 
easily  which  branch  of  the  multivalued  complex  function  one  must 
take. 

In  the  example  which  we  treat,  the  airfoil  profile  has  been 
scaled  so  that  the  trailing  edge  is  at  (-1,0)  and  the  point  (1,0) 
is  approximately  half  the  distance  between  the  leading  edge  and  the 
center  of  curvature  of  the  leading  edge.  The  portion  of  the  hor- 
izontal axis  between  -1  and  1 lies  inside  the  profile.  Figure  7. 
This  configuration  simplifies  the  mapping  procedure  somewhat. 

Highly  cambered  profiles  or  profiles  which  differ  considerably  from 


the  type  illustrated  may  require  a slightly  modified  treatment. 

Let  \ \ + iy^  for  k = 1 through  4.  The  Zj^-plane  is  the 

physical  plane,  that  is,  the  plane  of  the  airfoil  profile.  Let 

^2  “ (Bl) 

Then 

» 2/(1  - (B2) 

and 

lim  (1  + Zj^)/(1  - z ) = -1  (B3) 

z^-w 

Thus  the  image  of  the  point  at  infinity  in  the  z^-plane  under  the 
mapping  (Bl)  is  z^  = -1.  The  inverse  mapping  is 

^1  ° 

and 

dz^/dz^  = 2/(1  + z^)^  (B5) 

The  next  mapping  is  a rotation  which  causes  the  image  of  the 
upper  surface  of  the  profile  at  the  trailing  edge  to  be  tangent  to 
the  real  axis.  We  have 

z^  » z^  exp(-19^)  (B6) 

where  0^  is  the  angle  which  the  upper  surface  of  the  profile  makes 
with  the  Xj^-axls  at  the  trailing  edge.  Figure  7.  Then 
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dz./dz  = exp(-i6.) 

^ i.  A 

The  point  z^  = -1  goes  into  the  point 
z^  = expi(it  - 9^) 

The  inverse  mapping  is 

^2  ' ^3  exp(ie^) 

and 

dz2/dz^  = exp(ie^)  (BIO) 

Let  9g  denote  the  angle  which  the  lower  surface  of  the  profile 
makes  with  the  x^-axis,  Figure  7.  Our  object  is  to  eliminate  the 
sharp  corner  at  the  trailing  edge.  This  is  accomplished  by  the 
mapping 

Z4  = Z3®  (BID 

where  g = Tt/(9„  - 0.).  Then 
B A 

dzjdz^  * %z^  ^ (B12) 

The  point  z^  = expi(ir  - 0^)  goes  into  the  point 

z^^  = expin(n  - 6^)/(9g  " 6^)  (B13) 

The  inverse  mapping  is  given  by 

Z3  • *4^^*  (B14) 


(B7) 


(B8) 


(B9) 


57 


and 


dz3/dz^  = 

Finally  a linear  transformation 


(B15) 


(B16) 


Is  carried  out  which  maps  the  Image  of  the  point  at  infinity 

in  the  physical  plane,  into  infinity  in  the  z-plane  and  the  image 
of  the  trailing  edge,  z^  * 0,  Into  z = -1.  Then 

dz/dz^  - -2z^^/(z^  - z^^)^  (B17) 

The  inverse  transformation  is 


^4  ° ^44^^  ~ 


(B18) 


and 


dz^/dz  = -2z^^/(z  - 1)^ 

The  sequence  of  mappings  defined  by  equations  (Bl), 

and  (B16)  determine  z as  a function  of  z^.  Replacing  z^ 

value  in  terms  of  z^  and  so  on  gives 

(1  + z^)^  + (1  - Z3)®exp(lgTT) 

z * 

(1  + Z3)*  - (1  - Z3)*exp(lgii) 


(B19) 

(B6),  (BID 
by  its 


(B20) 


This  is  the  Karman  Trefftz-llke  transformation  z = FCz^)  which  we 
use  to  map  profiles  into  near  circles.  Figure  8 shows  the  near 
circle  image  of  the  Clark  Y profile  data  under  the  transformation 
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just  described. 


In  order  to  determine  the  values  of  U and  V,  Eq . (6)  the  value 
of  dF/dZj^  at  = “>  is  needed.  By  the  chain  rule 

dz/dz^  = (dz2/dz^) (dz^/dz2> (dz^/dz^) (dz/dz^)  (B21) 

As  z^  tends  to  infinity  dz/dz^^,  as  expressed  by  Eq  (B21),  is  inde- 
terminate, but  it  can  be  shovm  that  the  limiting  value  is  1/g.  Thus 
we  have  at  z^^  = ■» 

dF/dZj^  = 1/g  (B22) 

The  image  of  a highly  cambered  profile  under  a Karman  Trefftz 
transformation  may  still  be  close  to  a circle,  but  not  necessarllv 
close  to  a circle  with  center  at  the  origin  of  the  z-plane.  In  this 
event  it  is  practical  to  make  a judicious  shift  of  the  origin  to 
the  center  of  gravity  of  the  contour. 

The  Joukowski  transformation,  which  is  a special  case  of  the 
Karman  Trefftz  transformation,  applied  to  two  line  segments  shows. 
Figure  9,  the  effects  of  the  blow  up  of  one  profile  on  another. 

Figure  9,  panel  A,  shows  two  parallel  line  segments.  Panel  B shows 
the  image  of  these  line  segments  under  the  inverse  of  the  first 
Joukowski  transformation.  The  lower  line  segment  in  panel  A has 
been  unfolded  to  produce  the  unit  circle  in  panel  B.  Panel  C shows 
the  image  of  panel  B under  a transformation  of  the  form  Az  + B. 

The  endpoints  of  the  upper  contour  are  transformed  into  (-2,0)  and 
(2,0)  respectively  and  the  unit  circle  is  transformed  into  a circle. 


I 

i 
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Finally,  panel  D shows  the  Image  of  panel  C under  the  Inverse  of 
a second  Joukowskl  transformation.  This  time  the  upper  contour  of 
panel  C is  unfolded  into  a nearly  circular  contour  and  the  circle 
of  panel  C is  transformed  into  a near  circle  also. 

In  Figure  10  biplane  data,  given  in  reference  [11],  is  trans- 
formed into  two  near  circles  by  a composite  mapping  consisting  of 
two  Karman  Trefftz  transformations  and  three  linear  transformations. 
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APPENDIX  C 


NUMERICAL  DETERMINATION  OF  THE  DENSITIES 


The  unknovm  density  function  k(t)  or  the  tangential  velocity 
v(t)  is  determined  by  an  Integral  equation  of  the  form 


f(s)  = +wk(s)  + f k(t)K(s,  t)  dt  (Cl) 

C 

The  integrand  k(t)K(s,  t)  is  a periodic  function  of  t for  each  value 
of  s.  Such  Integrals  can  be  evaluated  with  high  accuracy  by  the 
trapezoidal  rule  provided  that  the  integrands  are  sufficiently 
smooth,  [ 7 ] . 

Since  the  contour  C is  described  as  r = T(6),  here  r denotes 
the  radius  vector  from  the  "center"  of  the  contour  C to  a point  on 
C,  9 instead  of  length  of  arc  t is  taken  for  the  integration  vari- 
able. Partition  the  interval  [0,  2it]  into  n equal  subintervals  and 
denote  the  endpoints  of  the  l*"^  subinterval  by  [0^, 
brevity  write  f^  for  f(9^),  for  K(6^,  0j)  and  so  on. 

Let  the  contour  C be  given  by  the  parametric  equations 

X = x(e) , y = y(0)  (C2) 


and  let  t(9)  denote  the  length  of  the  contour  C as  a function  of 
9.  Derivatives  of  functions  of  0 will  be  denoted  by  primes.  Then 
the  kernel 


K 


Ij 


(Xj  - - (y^^  - y^)x]^ 

(<x^  - Xj)^  + (y^  - 


(C3) 
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Is  Indeterminate  when  j “ i.  Applying  L' Hospital's  Rule  one  ob- 


tains for  the  limit  as  0j  tends  to  6^^ 


yVx* 


1-1  -17 1 


xVy; 


2((xp^  + (yp^]t* 


(C4) 


APPENDIX  D 


NUMERICAL  DETERMINATION  OF  THE  DERIVATIVES  OF  THE  DENSITIES 

The  derivatives  of  the  density  functions  are  needed  for  the 
numerical  determination  of  the  tangential  velocities,  Eq.  (b). 

These  derivatives  can  be  obtained  by  numerical  differentiation  of 
the  density  functions.  They  can  be  obtained  also  from  the  Integral 
equation  as  shown  in  this  appendix. 

Differentiating  Eq.  (Cl)  with  respect  to  arc  length  s gives 

df/ds  = TT(dk/ds)  + f k(t)K„(s,  t)  dt  (Dl) 

C 

where  “ 3K(s,  t)/3s.  Set 

A (x  - Xj)(d^y/ds^)  - (y  - y^)  (d^x/ds^) , 

B = (x  - Xj)(dy/ds)  - (y  - yj)(dx/ds), 

C = (x  - Xj)(dx/ds)  + (y  - yj)(dy/ds), 

R = (x  - 4 (y  - y^)^ 

Then 

K2(s,  t)  = (RA  - 2BC)/R^ 

The  kernel  K2(s,  t)  is  indeterminate  for  t = s. 

A 4 4 

Let  us  write,  for  example,  R for  d R/dt  . Using  this  conven- 
tion, one  has  after  four  applications  of  L'Hospital's  Rule 

4 31  22  13  4 4 31  22  13  4 

„ , , RA  4 4RA  4 6RA  4 4RA  4 RA  - 2(BC  4 4BC  4 6BC  4 4BC  4 BC) 

K2(s.s) 22 13 4 

6RR  4 8RR  4 2RR 
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which  simplifies  to 


13  31  22 

K2(s,  s)  = (1/3) (RA  - 2BC  - 3RC)  (D2) 

The  right-hand  side  of  Eq . (D2)  written  out  is 

(2/3) ( (dx/ds) (d^y/ds^)  - (dy/ds) (d^x/ds^)) + 

((dx/ds) (d^x/ds^)  + (dy/ds)(d^y/ds^))((dy/ds)(d^x/ds^)  - (D3) 

(dx/ds) (d^y/ds^)) 


We  are  now  ready  to  evaluate  numerically  the  integral  in  Eq . (Dl) . 
We  abide  by  the  same  conventions  as  in  Appendix  C and  take  0 as  the 
Integration  variable.  Then 


n 

kjK2j^(dt/d0)j 

for  i = 1,  ...  , n.  Now  (dk/ds)^  is  the  value  of  dk/ds  at 
tlon  0^.  Solving  Eq . (D4)  for  (dk/ds)^  and  multiplying  by 
gives  the  value  of  dk/d0  at  the  station  0^. 


(df/ds)  = ir(dk/ds)  + A0  X 

j-1 


(D4) 


the  sta- 
(ds/d0)^ 
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APPENDIX  E 


NUMERICAL  DETERMINATION  OF  THE  TANGENTIAL  VELOCITIES 


In  this  appendix  we  discuss  the  evaluation  of  the  integral  in 
Eq . (5)  for  the  tangential  velocities.  The  equation  for  the  tan- 
gential velocities  is  quite  similar  to  Eqs.  (Cl)  and  (Dl).  The 
difference  is  that  the  kernel  in  Eq . (5)  has  a first  order  singu- 
larity for  t = s.  Hence  the  integral  must  be  interpreted  in  the 
sense  of  its  Cauchy  principle  value. 

Denote  the  integrand  in  Eq.  (5)  by  H(6).  We  have 

H(e)  = (x(e)-x(9))/(dx(e)/dt)-f(y(9)-y(9))(dy(e)/dt)  j,(e)(dt/d9) 
(x(9)  - x(9))^  + (y(9)  - y(9))^ 


and  H(9)  has  a first  order  singularity  at  9 = 9.  Take  't'  •=  9 - 9. 
Then  we  can  write 


H(9)  = (1/2)  [H(9+((i)  + H(9-4i)]  + (1/2)  [H(S+<ti)  - H(9-i(.)]  (El) 


that  is,  as  the  sum  of  an  even  and  odd  function  of  <>.  It  is  clear 

that  the  singularity  is  in  the  odd  part  of  H(9).  Then 

2n  n , 

CPVy"  H(9)d9  = J (l/2)[H(9+(ti)  +H(9-()i)]  d<)>  (E2) 

^0  -TI 

since  the  Cauchy  principle  value  of  the  integral  of  the  odd  part  of 
H(9)  is  zero. 

Needed  for  the  evaluation  of  the  right-hand  side  of  Eq.  (E2)  is 
the  value  of 


lim  (1/2)[H(9  + 1^)  + H(9  - 4>) ) 


If  H(0)  were  a continuous  function  of  9 then,  as  noted  above,  we 
would  have 


2tt 

f H(0)  d0  = A0  Z H(0  ) (EA) 

0 j J 

Suppose  now  that  the  first  order  singularity  of  H(9)  Is  at  the  i*'*^ 
station  0^.  Then  the  sum 

AOZhOJ  + A0  lim  (1/2)[H(9  + e)  + H(0.  - e)]  (E5) 

j J e-0  ^ ^ 

where  Z'  denotes  the  sum  over  all  values  of  j i,  is  equivalent 

to  evaluating  the  even  part,  that  is,  the  right-hand  side  of  Eq. 

(E2),  by  the  trapezoidal  rule. 

Set 

A = (dx/dt) (dx/d0)  + (dy/dt) (dy/d0)  (E6) 

B = k(9) (d^t/d9^)  + (dk/d0)(dt/d0)  (E7) 

C = 0.5( (dx/dt) (d^x/d0^)  + (dy/dt) (d^y/d0^)]k(9)(dt/d6)  (E8) 

D = (dx/d9)^  + (dy/d0)^  (E9) 

E - ( (dx/d0)(d^x/d0^)  + (dy/d0)(d^y/d0^)]k(0)(dt/d0)  (ElO) 

Then  the  limiting  value  of  O.5[H(0  + <>)  + H(0  - i>)]  as  tends  to 
zero  is  given  by  the  expression 

-(AB  + C)/D  + A£/D^  (Ell) 

evaluated  at  0. 
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APPENDIX  F 


CURVE  FITTING  AND  SMOOTHING 


We  have  assumed  that  the  contour  C is  given  as  a function  of 
a continuous  parameter.  In  practice,  the  airfoil  profile  frequently 
is  given  as  a discrete  set  of  points.  In  this  case  some  curve  fit- 
ting process  is  needed  to  determine  the  profile  as  a function  of  a 
continuous  variable.  After  the  profile  data  is  transformed  to  the 
near  circle  plane  the  rapid  change  at  the  leading  edge  and  the 
sudden  change  at  the  trailing  edge  have  been  removed.  For  these 
reasons  it  is  expected  that  the  near  circle  data  can  be  interpolated 
easier  than  the  profile  data.  We  use  a cubic  spline  to  fit  the  near 
circle  data.  The  general  equation  for  the  cubic  spline  and  the  pro- 
cedure for  fitting  the  near  circle  data  is  given  next. 

Let  there  be  given  n data  points  (x^^.y^)  with  x^  < We 

approximate  the  curve  between  two  points  (x^,y^)  and 
with  a polynomial  p^(x)  of  the  third  degree  and  of  the  form 


^ ^"i+1  ^1  v^i  /j+i 

° 6(xi^l  - x^)  “ 6(x^^^  - x^) 

+ - x)  - x^)/6)  (Fl) 

+ (x  - - x^)  - yj'+i^^i+i  ■ 


where  the  coefficients  are  expressed  in  terms  of  the  quantities  y'^ 

and  y'^^2^"  The  resulting  curve  has  continuous  curvature  over  the 

whole  Interval  (x,  , x ) . 

1’  n 
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In  the  present  case  y Is  a periodic  function  of  x.  Using 


this  fact  we  obtain  the  following  system  of  equations  for  the 


s : 


(x 


1+1 


xi)yi’  + 


2(x 


1+2 


- ’'i^yi+i 


+ (x 


1+2 


1+1 


>^1+2 


^(^1+2  - ^1+1  _ yj+1  - yj 

’'i+2  " ’'i+1  ^1+1  ■ ’'l 


(F2) 


for  1-1,  ...  , n-3.  For  1 = n - 2 and  n - 1 we  have,  respectively 


^^n  - ’‘n-l>yi  + ^’‘n-l  " ’'n-2>yn-2 
(.(h-Llnzl  _ ^n-l  ~ yn-2 

*n  ■ ’‘n-l  ’‘n-1  ~ ^n-2 


(F3) 


2 (x  - X , + 
n n-1 


^2 


xi)y'i'  + 


(X2  - 


xi)y^’  + 


(x  - 
n 


X ,)y" 
n-1  ^n- 


yi 


X 

n 


•) 


(F4) 


Note  that  because  of  the  periodicity  we  have  n-1  Instead  of  n 
equations  and  unknowns. 

The  airfoil  profile  Is  transformed  into  a near  circle  in  the 

z-plane.  If  the  profile  is  given  by  a set  of  n discrete  points  then 

the  near  circle  consists  of  a set  of  discrete  points  z^,  1 ••  1,  ...  , 

n with  z » z, . Let  r . and  0,  denote  the  polar  coordinates  of  the 
n 1 11 

point  z^.  Taking  6 as  the  Independent  variable  we  set  up  the  sys- 
tem of  equations  (F2),  (F3)  and  (FA)  with  r^  playing  the  role  of 
y^  and  0^^  playing  the  role  of  x^. 
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In  Figure  3 we  show  the  cubic  spline  fit  of  the  near  circle 
data  transformed  into  the  physical  plane  by  the  Inverse  of  the  Kar- 
man  Trefftz-like  transformation.  First  impressions,  from  this  fig- 
ure, are  that  the  cubic  spline  fit  is  more  than  adequate.  However, 
the  graph  of  the  derivative  of  this  cubic  spline  fit.  Figure  4, 
shows  considerable  wavlness.  Thus  inaccuracies  in  the  profile  data 
may  produce  rapidly  varying  slopes  and,  in  a high  precision  method, 
the  pressure  results  will  reflect  these  variations  in  slope.  This 
may  be  the  reason  practitioners  recommend  limited  accuracy  in  such 
computations.  Obviously,  the  wavlness  of  the  results  obtained  with 
high  accuracy  computations  is  not  caused  by  the  method  but  by  rough- 
ness in  the  data  to  which  the  method  is  applied.  We  mention  this 
fact,  because  if  it  is  not  clearly  recognized  a good  method  might 
be  rejected  because  the  results  are  sensitive  to  inadequacies  in  the 
data . 

Smoothing  the  pressure  results  caused  by  inaccuracies  in  the 
profile  data  Introduces  an  element  of  arbitrariness.  It  seems  pre- 
ferable to  smooth  the  profile  data.  One  can  proceed  as  follows: 
After  r is  determined  as  a function  of  the  continuous  variable  9, 
the  values  of  dr/d0  are  computed  and  plotted.  If  the  graph  of  dr/d6 
oscillates  too  much  some  of  these  oscillations  should  be  smoothed 
out . 

The  graphs  of  r - 1 versus  9 and  dr/d9  versus  9 are  given  in 
Figure  4.  Using  a commercial  reader  (a  digitizer)  values  of  dr/d9 
are  read  at  equal  increments  of  9 from  a smoothed  graph  of  dr/d9. 
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These  values  of  dr/d6  are  used  to  determine  a trigonometric  inter- 
polating polynomial  for  dr/d0.  From  this  trigonometric  interpolating 
polynomial  r is  redetermined  as  a function  of  0,  r * T(0).  Figure 
5 shows  the  gr.iph  of  dr/d0  (not  smoothed),  the  values  of  dr/d0  which 
were  read  at  equal  Increments  of  0 from  a smoo'.hed  graph  of  dr/d0  and 
the  smoothed  graph  of  dr/d0  as  given  by  the  trigonometric  inter- 
polating polynomial.  The  smoothed  graph  as  given  by  the  trigono- 
metric interpolating  polynomial  does  not  pass  through  the  smoothed 
data  because  the  constant  term  must  be  set  equal  to  zero.  (Other- 
wise the  curve  will  not  close  to  the  same  value  of  r after  one  com- 
plete cycle.) 

To  check  the  smoothing  process  the  points  (T(0),0)  for  0^0^ 
2tt  are  transformed  to  the  physical  plane  by  the  inverse  of  the  Kar- 
man  Trefftz  transformation.  In  this  way  a new  airfoil  profile  is 
determined.  If  this  new  profile  is  technically  unacceptable  then 
the  smoothed  values  of  dr/d0  need  adjusting.  In  Figure  2 we  show 
the  graph  of  the  image  of  (T(0),0)  for  0 ^ 0 ^ 2ir  in  the  physical 
plane  and  the  plot  of  the  original  data. 

The  pressure  distribution  is  computed  for  this  new  profile. 

Of  course,  if  the  fitting  of  the  original  data  and  the  smoothing  is 
performed  with  sufficient  skill,  the  new  profile  will  be  close  to 
the  original  profile.  After  the  fitting  and  smoothing  process  is 
satisfactory  the  next  step  is  the  determination  of  the  density 


APPENDIX  G 


THE  COMPUTER  PROGRAMS  AND  LISTINGS 

In  the  first  part  of  this  appendix  we  review  the  conditions 
which  must  be  satisfied  in  order  to  use  the  computer  programs  listed 
at  the  end.  These  conditions  and  Instructions  are  given  in  the 
comment  statements  also,  SECTIONS  II,  III  and  the  preceding  appen- 
dices provide  the  theoretical  basis  and  details  for  the  creation 
of  these  programs.  These  programs  fall  into  two  categories.  In 
the  first  category  are  four  programs  for  determining  contributions 
to  the  tangential  velocities  (usually  in  the  plane  of  the  near  cir- 
cles). In  the  second  are  programs  concerned  with  mapping  an  air- 
foil profile  into  a near  circle. 

The  program  which  maps  an  airfoil  profile  into  a near  circle 
is  the  first  step  in  determining  the  equations  for  the  near  circle 
contour.  The  equations  for  the  near  circle  contour  (or  contours) 
are  needed  in  the  programs  for  computing  the  contributions  to  the 
tangential  velocities.  The  point  which  we  wish  to  make  here  is  that 
the  listed  programs  only  provide  intermediate  results.  These  are 
working  programs  but  they  are  not  finished  products. 

The  four  programs  MTVLPI,  MTVELP,  MTVLSI  and  MTVELS  for  the 
contributions  to  the  tangential  velocities  are  much  alike  and  in- 
structions for  using  them  are  the  same  in  each  case.  The  programs 
MTVLPI  and  MTVELP  are  based  on  Prager's  method  with  vortices  dis- 
tributed along  the  contours.  The  programs  MTVLSI  and  MTVELS  assume 
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a distribution  of  sources  along  the  contours.  The  programs  MTVLPl 
and  MTVLSI  solve  certain  systems  of  equations  iteratively.  Note 
that  the  letter  S in  the  program  name  Indicates  a distribution  of 
sources  around  the  contours,  the  letter  P a distribution  of  vortices 
and  the  letter  1 means  that  systems  of  linear  equations  are  solved 
by  an  iterative  process. 

These  programs  were  developed  on  the  CDC  Cyber  74.  The  pro- 
gramming is  standard  FORTRAN  IV  except  for  the  use  of  NAMELIST  in 
conjunction  with  the  RE7vD  and  WRITE  statements  on  unit  5.  NAMELIST 
is  used  only  for  separating  and  locating  records  on  the  file  asso- 
ciated with  unit  5.  These  nonstandard  statements  occur  only  in  the 
main  or  control  programs,  that  is,  in  PROGRAMS  MTVLPl,  MTVELP,  MTVLSI 
and  MTVELS.  A user  will  have  to  replace  the  programming  noted  with 
programming  which  is  compatible  with  the  available  processor  and 
which  provides  the  capability  of  reading  previously  prepared  records 
and  writing  new  records.  Also  in  developing  these  programs  it  was 
convenient  to  write  formatted  records.  For  general  use  format  free 
records  are  probably  preferable. 

The  information  which  must  be  supplied  to  these  programs  is 
the  number  of  contours  NC,  the  number  of  points  N(I)  around  the  1th 
contour  at  which  contributions  to  the  tangential  velocities  are  to 
be  computed  and  the  equations  of  the  near  circle  contours.  The 
equations  of  the  near  circle  contours  are  given  in  terms  of  a local 
polar  coordinate  system.  The  radius  vector  from  the  (approximate) 
center  of  the  near  circle  contour  to  a point  on  the  contour  is 
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assumed  given  by  a trigonometric  interpolating  polynomial.  Thus 
the  near  circle  contours  are  determined  by  their  (approximate)  cen- 
ters and  the  coefficients  of  the  trigonometric  interpolating  poly- 
nomials. 

The  integers  NC  and  N(I)  for  1=1,  ...  , NC  are  supplied  at 
input,  F0RMAT(I5).  Using  NAMELIST  to  locate  the  proper  record,  the 
coefficient  number  M(J),  F0RMAT(I5),  the  coordinates  (XX(J),  YY(J)) 
of  the  "center"  and  the  coefficients  A(I,  J) , for  1=1,  ...  , 

M(J)+1,  B(I,  J)  for  1=1,  ...  , M(J)-1  of  the  trigonometric  inter- 
polating polynomial  of  the  Jth  contour,  FORMAT(5E23.14) , for  J = 1, 
...  , NC  are  read  from  the  file  associated  with  unit  5.  A user  may 
have  to  provide  some  alternate  way  to  supply  the  contour  equations 
to  these  programs. 

The  principle  outputs  from  the  programs  MTVLPI  and  MTVELP  are 

the  arrays  AA(NT),  BB(NT)  and  CC(NT,NC).  Here  NT  =J]n(I).  The 

I 

tangential  velocity  at  points  on  the  near  circle  contours  is  a linear 
combination  of  the  entries  in  the  arrays  AA,  BB,  and  CC . The  co- 
efficients for  Liie  linear  combination  are  determined  by  the  angle 
of  attack  and  the  Kutta  condition.  In  the  programs  MTVLSI  and 
MTVELS  the  tangential  velocity  is  obtained  as  a linear  combination 
of  the  entries  in  the  arrays  TAA(NT) , TBB(NT)  and  TCC(NT,  NC) . 

These  arrays  are  written  on  the  file  associated  with  unit  5.  The 
user  may  have  to  provide  another  way  for  saving  and  recalling 
this  information  for  computing  tangential  velocities  and  the  pres- 
sure coefficient. 
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The  program  for  mapping  profile  data  into  near  circle  data  by 
a Karman  Trefftz-like  transformation  is  called  KATR.  The  Karman 
Trefftz-llke  transformation  is  multivalued.  The  program  listing 
of  KATR  is  Intended  only  for  profiles  which  have  the  form  indicated 
by  Figure  7.  That  is,  if  one  draws  a horizontal  line  through  the 
trailing  edge,  then,  except  for  the  leading  and  trailing  edges,  the 
upper  surface  of  the  profile  lies  entirely  above  this  horizontal 
line  and  the  lower  surface  lies  entirely  below.  The  Karman  Trefftz 
transformation  is  not  limited  to  profiles  for  which  the  upper  and 
lower  surfaces  satisfy  the  above  condition. 

The  profile  data  is  separated  into  two  sets.  One  set  consists 
of  the  points  on  the  upper  and  the  other,  the  points  on  the  lower 
surface.  The  leading  and  trailing  edges  are  Included  in  both  sets. 
NA(NB)  denotes  the  number  of  data  points  on  the  upper (lower)  sur- 
face. Beginning  with  the  leading  edge,  the  points  on  the  upper 
(lower)  surface  are  indexed  from  1 to  NA(NB) . It  is  assumed  that 
the  profile  data  is  scaled  so  that  the  trailing  edge  has  coordinates 
(-1,  0).  The  distance  from  the  point  (1,0)  to  the  leading  edge  is 
approximately  half  the  radius  of  curvature  at  the  leading  edge. 

The  input  data  to  the  program  KATR  is  the  Integers  NA  and  NB, 
FORMAT(215).  Next  the  coordinates,  FORMAT ( 6E10. 0) , of  the  points 
on  the  upper  and  lower  surfaces  are  read  from  input.  We  used  NAMELIST 
in  conjunction  with  WRITE  statements  (as  noted  above)  to  place  this 
input  information  on  the  permanent  file  associated  with  unit  5.  A 
user  may  have  to  replace  this  coding  with  coding  compatible  with 
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his  central  processor. 

The  output  from  KATR  is  the  four  parameters  TA(=0.),  TB(=6  ),G 

A O 

and  T which  are  used  in  defining  the  (Carman  Tref f tz-like  transfor- 
mation and  its  inverse,  and  the  entries  in  the  arrays  TH(K),  RO(K) 
and  R2(K),  where  K = NA+NB-1.  The  entries  in  the  arrays  RO  and  TH 
are,  respectively,  the  polar  coordinates,  radius  vector  and  polar 
angle,  of  the  near  circle  images  of  the  profile  data.  The  entries 
in  the  array  R2  are  the  second  derivatives  of  the  cubic  spline  fit 
(the  second  derivatives  of  the  radius  vector  with  respect  to  the 
polar  angle)  of  the  near  circle  data.  RO(K)  = R0(1),  R2(K)  = R2(l) 
and  TH(K)  » TH(1)  + 2tt. 
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PROG-<AM  MTVLPI  ( T APE  5 , 1 NPU  T , OU  TPUT  ) 

01  HE  NS  I ON  AAA1(2,2),IP1(2),AA1(2),B31(2) 

CONN ON/ JK3/i 1( 3 J) , 32(90), X(0O),X1(3Q), 

2 X2(90)  ,Y(30),Yl(a0)  ,Y2(a0),XlP(80), 

3 YIP (80 ) ,T I ( 2) ,SO(80) 

CONNON/ 3K4/AA(90),33(80),CC(90,2),CI(80),CA(80,2) 
CONNON/ 3K5/AAA (90,90),IP(80) 

CONNON/ JKo/N(2) ,NC,NT 

CONNON/  3K7/NK2)  ,XNE  (2)  ,NA(2)  ,XAM(2) 

CONNON/3K3/XO(2),YO(2),M(2),A(51,2),0(49,2) 

♦*•♦*♦•*•♦•***♦»♦*•****•*♦♦*♦•♦♦***»•♦♦*♦*♦*»♦* 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


NTVLPI  IS  THE  CONTROL  PROGRAM  FOR  DETERMIN- 
ING THc  tangential  \/EL0CITIES  INDUCEJ  3Y  A 
0ISTRI3JTI0N  OF  VORTICES  ON  THE  CONTOURS. 

IN  this  program  most  of  THE  LINEAR  SYSTEMS 
OF  EQUATIONS  ARE  SOLVED  3Y  ITERATIVE 
PROCESSES. 

NC  NUM3ER  OF  CONTOURS 

N(I)  NUN3ER  OF  SU3DIVISI0NS  OF  THE  ITH  CONTOUR 

NT  total  number  OF  SUBDIVISIONS  ALl  CONTOURS 

NT=SUM  N ( I) , 1=1, NC 

WE  HAVE  USED  PIXED  RATHER  THAN  VARIABLE  DIMENSIONS.  IF  THE 
NUMBER  OF  CONTOURS  OR  IF  THE  NUMBER  OF  SUBDIVISIONS  IS  I N- 
CREASED  Then  the  DIMENSION  STATEMENTS  WILL  HAVE  TO  BE  CHANGED 
AND  The  program  MTVLPI  AND  SUBROUTINtS  WHICH  IT  CALLS  WILL 
HAVE  TO  3E  RECOMPILED  TO  REFLECT  THE  NEW  VALUES  OF  NT  AND  NC 
XO(K),YO(K)  APPROXIMATE  "CENTER"  OF  THE  KTH  CONTOUR 
A (L ,K)  ,0 (L  ,K)  FOURIER  COSINE  AND  SINE  COEFFICIENTS  FOR  THE 

RADIUS  VECTOR  R FROM  XO(K),YO(K)  TO  A POINT  ON 
THE  KTH  contour. 

the  centers  ( X0(<) , Y0(K) ) AND  THE  FOURIER 
COEFFICIENTS  A(L,K)  AND  3(L,K)  ARE  READ 
FROM  the  file,  TAPE5.  THESE  QUANTITIES 
MUST  3i  PREPARiD  ANO  STORED  IN  ADVANCE 
IN  THE  PROPER  FORMAT. 

AA,BB  CONTRIBUTIONS  TO  THE  TANGtNTIAL  VELOCITY 

R-SULTING  FROM  THE  VELOCITY  COMPONENTS 
AT  INFINITY. 

CC  COMPONENTS  OF  THE  EIGENVECTORS.  THESE 

PROVIDE  the  CIRCULATION. 

CA  THE  EIGENVECTORS  OF  THE  ADJOINT 

PROBLEM  OR  FOR  THE  TRANSPOSED  MATRIX 
X(I),Y(I»  COORDINATES  OF  THE  ITh  STATION  ON  THE  CON- 

TOURS. XKDjYKI)  ARE  THE  FIRST  DERIVATIVES 

w.r.t.the  central  angle  theta  and  so  on 

XIP(I) ,Y1P(I)  ARE  THE  FIRST  DERIVATIVES 
W.R.T.  arc  length  and  SO  ON 
S1,S2  denote  first  ANO  SECOND  DERIVATIVES 

OF  THE  ARC  length  W.R.T.  THE  CENTRAL  ANGLE 

SO  Change  of  integration  variable  factor 

MTVLPI  CALLS  THE  FOLLOWING  SUBPROGRAMS 

MVALP  THIS  SUBROUTINE  USES  THE  QUANTITIES  IN  COMMON  BLOCK 
a TO  COMPUTE  The  quantities, PRIMARILY  IN  BK3. 
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c 

c 

c 

c 

c 

c 

r. 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


'^P'iAri  THIS  SJBr^OUri'JE  USES  THE  VALUES  IN  COMMON  BLOCK  J 
TO  COMPUTE  The  ENTRIES  FOR  THE  MATRIX  AAA,  THt-N 
THt.  EIGENVECTORS  GORrESPONO  ING  TO  THE  EIGEN- 
VALUE PI  OF  AAA  AND  ITS  TRANSPOSE  ARE  COMPUTED 
THc:  MATRIX  ENTRIES  ARE  STORED  IN  COMMON  BLOCK  5, 
SCV  THIS  subroutine  GIVES  THE  PRODUCT  OF  A VECTOR 

BY  A scalar. 

PMAT2  THIS  SUBROUTINE  COMPUTES  THE  COEFFICIENT 
MATRIX  USlO  in  REMOVING  THE  COMPONENTS  IN 
IN  THE  OiRECTION  OF  THE  EIGENVECTORS  FROM 
A GIVEN  VECTOR. 

LUO  A ROJTINl  for  THE  DIRECT  SOLUTION  OF  SYSTEMS 

OF  LINEAR  EQUATIONS  TAKEN  FROM  THE  ARL 
LINEAR  algebra  LIBRARY.  IT  INCLUDES  LEQS3. 

PITER  this  subroutine  SOLVlS  FOR  THE  CONTRIBUTIONS  TO 
the  tangential  VELOCITIY  by  an  ITERATION 
PROCESS.  THIS  ROUTINE  IS  APPLICABLE  ONLY 
TO  VERY  SPECIAL  SYSTEMS  OF  LINEAR  EQUATIONS 
the  input  data  TO  RUN  THIS  PROGRA M IS  AS  FOLLOWS 
NOTE-THE  LETTER  d DENOTES  A BLANK  IN  THE  INPUT  OATA 
♦EOR 
B33NC 
BBBNd  I 
ETC. 

OBBN(NC) 

*EOR 


NAMELIST/NAM/LLL/NAM  l/LLL/NAM2/LLL/NAM3/LLL/NAM‘f/LLL/NAM5/ 
2 LLL^UAMb/LLLYNAH7/LLL/NAMft/LLL/NAMD/LLL 
LLL  = 1 

PI  = 3. 14  15D25535398 
REAj  B00,NC 
NT  = 0 

□ 0 10  1=1, NO 
read  9Q0,N(I) 

10  NT  = NT  + Na> 

RE  AO (B,  NAM4) 

DO  20  J=1,NC 
READ (5, BOO ) M ( J) 

MPr.M  ( J)  <-1 
HM=M ( J) -1 

READ  (5  , TOl ) XO( J)  , YO( U> 

REA0(5,  301) (A(I,J),I=1,MP) 

20  REAJ(5,  301)  (t3(  I ,D)  , r = l,MM) 
call  MVALP 
C=1/PI 

PI  call  SCV(NT,C,AA) 

CALL  SCV(NT,C,BB) 
call  MOMATI 

t.  any  :J3P0NENTS  along  the  EIGlNV^CTORS  CORRESPONDING 
f ) T.-*E  lIGlNVALUl  °l  ARE  REMOVED  NEXT. 

CAl.  PMATe (NT, NC,AAA1,CA,CC) 

0)  ..’4  jd,N: 

AA I ( n = ) 


78 


asi(j»=Q 

DO  22  1=1, NT 

Aftl(J)=AAl(J)fCA(I,J)*AA{I) 

22  8B1(J)= 3B1(J)+CA(I, J)*3B(I) 

2^  CONTINUE 

call  LU3(NC,AAA1,2,IP1) 

CALL  LE 3S3 (NC, AAAI, 2 , IPl, AAl) 
call  LEUS3 (NC, AAA1,2 ,IP1,3B1) 

DO  2b  1=1, NT 
DO  25  J=1,NC 

AA(I)=AA(I)-AA1(J)*CC(I,J) 

25  B3(I>  =8-J<r  ) -031  ( J)  *CC  (I , J) 

26  continue 

C THE  COMPONENTS  IF  ANY  ARE  REMOYEJ 
WRITE (5  ,NAM5 ) 

C STORE  THl  EI3ENVECTOR5  OF  THE  TRANSPOSE  OF  AAA. 
DO  40  <1=1, NC 

11  = 1 
12=0 

WRI TE (5 ,300) NA ( Kl) 

WRITE (5,901) XAH (<1) 

00  30  K=1,NC 
I2  = I2  + N (<) 

WRITE (5 ,900) N<  K) 

WRITE(5,90l) (CA(I,<1),I=I1,I2) 

30  I1=I1^N(K) 

40  continue 

WRITE (5  ,NAHb) 

call  PITER(NT, K,EH,AAA,AA,L) 

IF  (L)  42,42,44 

42  PRINT  902, EH 

GO  TO  100 

44  WRITE (5 ,900) K 

WRITE (5  ,90 1) EH 
11=1 
12=0 

DO  50  K=1,NC 
I2  = I2  + N (K) 

WRITE(5  ,900) N(<) 

WRITE(5,901)  (A  A ( I) , I =I 1, 12) 

50  H = I1  + N(K) 

call  PITER(NT,K,EH, AAA,3B,L) 

IF  (L)  52,52,54 

52  PRINT  902, EH 

GO  TO  100 

54  WRITE (5  ,900) < 

WRITE(5  ,901) EH 
11=1 
12  = 0 

00  oO  K=1,NC 
I2  = I2  + N (K) 

WRITE(5,900) N(<) 

WRITE (5, 90  1)  (B3(I),I  = I1,I2) 

60  I1=I14N(K) 
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03  ia  Kl=l,NC 
11  = 1 
1^=3 

WRITE  13  ,'iJU)  NI  (<1) 

WRITE.  (?  ,011)  XWE  (Kl) 

03  i^O  K=1,NC 
I2  = I2  + N (K> 

WRITE (5 ,930) N(K) 

WRIT£(5,90  1)  (CC(I,K1), 1 = 11,12) 

70  I1=-1+N(K) 

80  CONTINUE 

100  STOP 

900  E0R9ftT(I5) 

901  F0R9AT  ( 3E23.  1*+) 

902  FOR'IAT  ( *OEM=*- 23,  14,  ♦,  ITERATION  PROCESS  010  NOT  CONVERGE*) 
EN3 


SJIRO'JT  INE  "IVAlo 
01 'TENS  I ON  CH(2) 

COIN  ON/  OKS/SKaO)  ,32(80)  ,X(80)  ,X1(80)  , 

2 X2(B0)  ,Y(30),Yl(80),r2(30),XlP(80), 

3 Yl^>(80  ),  ri(  2)  ,SD(90) 

CJ990N/  ■3K4/AA(80),  33  (b0),CC(90,2),CI(80),CA(63,2) 
COMMON/ 3K8/N (2 ) ,NC,NT 

C09M0N/ 3<8/XO( 2 ) , YO ( 2)  ,R ( 2)  ,A (51 ,2) , b(49,2) 
real  NFOU, 93F0U ,M02F0U 


c 

c 

c 

c 

c 

c 

c 

r 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


this  SJ3ROUTINC  PREPARES  TABLES  OF  VALUEb  TO  BE  USED 
IN  LATER  COMMUTATIONS.  SPECIFICALLY,  IT  COMPUTES  THE 
C3OR0INATES  AT  STATIONS  ALONG  THE  CONTOUR  C,THE 
DERIVATIVES  OF  THESE  COORDINATES  W.R.T.  THE 
POLAR  variable  THETA=T,AN3  WITH  RESPECT  TO  ARC 
length.  these  AR-  denoted  E.G.,BY  XI  ANO  X1P,RESP. 

IT  Al33  computes  THE  DERIVATIVES  OF  T 3E  ARC  LENGTH 
W.R.T.  THE  POLAR  VARIABLE.  THESE  DERIVATIVES  ARE 
DENOTEj  by  S1,S2,ETC. 

This  subroutine  uses  thc  following  real 
function  PROGRAMS! 

MFOU  EVALUATES  A TRIGONOMETRIC  POLYNOMIAL  WITH 
FOURIER  COEFFIC IdNTS  A(I,K)  AND  IUI,K)  AT 
A olVEN  angle  (IN  RADIANS). 

MDFOU  EVALUATES  THE  FIRST  DERIVATIVE  OF  THE 
T R IGONOME TRIC  POLYNOMIAL. 

MD2F0U  EVALUATES  THE  SECOND  DERIVATIVE  OF  THE 
TRIGONOMETRIC  POLYNOMIAL. 


PI=i. 1415926535898 
DO  10  K=1,NC 
10  CH(O=0. 

11  = 0 

DO  40  K=l, NC 
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TI (<) =2 .*PI/N( K) 

T=-TI (K) 

I2=N(K) 

00  JO  1=1,12 
T = TfTI  ( <) 

R = iF0U(T,'1(<),/V,3,<) 
Rl=iOFO'J(T,M  (K)  ,A,i3,K) 

R2  = ^i02FaU(T,M(<)  ,A,3,K) 
Sl(I«-II)  = SaRT(R**2  + Rl**2) 
SOII+II  )=S1  ( I4-II)  *TI  (K) 

U = R«-R2 

y/=Ri*u 

S2(I»II )=7/Sl( If  II) 

02=S1 (I fll ) f*2 
C0=G0S ( T) 

SI=SIN( T) 

X(Ifll) =COfRfXO(K) 

Y(Ifll) =SI*RfYO(K) 

Xl(IfIl)=C:0*Rl-RfSI 

YKIfll  )=3I*RlfRfC0 

X2(IfII )=G0*R2-2. *31 *Rl-C0fR 

Y2 ( If  II )=SI*R2f 2.*C0*R1-SI*R 

XlP(IfII)=Xl(IfII)/Sl(IfII) 

Aft(ifII)=2.*PI*XlP(IfII) 

YlP(IfII)=Yl(IfII)/Sl(IfII) 

BB(IflI)=2.*PI*YlP(IfII) 

30  C0MTINU£ 

Il  = IIfN  (K) 

40  CONTINUE 

DO  70  K=1,NC 
11=1 
12=0 
CR(<) =1 
DO  oO  J=1,NC 
I2=I2fN  (J> 

00  50  1=11,12 
50  CC(I,K) =CH(J) 

60  Il=I2fl 

70  CH(O=0 

return 
END 


SUOROUTINi.  SCV  (N,C,  V» 
dimension  V(80) 


C this  SU3ROJTINC  MULTIPLIES  A VECTOR  BY  A SCALAR. 

00  10  1=1, NT 
10  V(I)=C*V<I) 

RETURN 

END 
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SU3^0UTIN-_  3P3ATI 

CaM-l3N/}K3/Sl(a0),S2  (80),X(60),Xl<30), 

2 X2(5C) ,Y(3O),yi(03),Y2(aQ),XlP(aO), 

J Y1P(80  ) , r I ( 2)  , 30( 30 ) 

C0  3 3 0N/jK-./aa(30»,3in8J),CC(30,2),CI(80),CA(80,2) 

C0M30N/ 3K5/flAA (30,30) ,IP(eO) 

COMMON/ 3K3/^ (2) ,NC,NT 

COMMON/  3K7/rJI(2),X3E;(2),NA(2),XAM(2) 

**««**»«»**** 

C THIS  3'J3ROiJTINl  JETERHIHES  THE  MATRIX  ELEMENTS  NEEOEO  FOR  THE 
C COMPUTATION  OF  THE  TANGENTIAL  VELOCITIES. 

C THEN  THE  EIGENFUNCTIONS  ARE  OETERMINED, 


I2=NT-1 
DO  30  1=1,12 
Jl=If 1 
DO  20 

U=(X(I)-X(J))*Y1P(I)-(Y(I)-Y(J))*X1P(I) 
V=(X(I)-X(J))**2«-(r(I)-Y(J))**2 
20  AAA ( I , J )= ( U/ V) *SD ( J) 

30  CONTINUE 

DO  30  1=2, NT 
J2=r-i 

00  >.0  U =1,  J2 

U=(X(I)-X(J))*Y1P(I)-(Y(I)-Y(J))*X1P(I) 
V=(X(I)-X(J))»*2+-(Y(I)-Y(J))**2 
40  AAA( I , J > = (U/ V)  *SO ( J) 

50  CONTINUE 

00  50  I =1, NT 

U=Y2(I)*X1P(I)-X2(I)*Y1P(1) 

V = X1  ( I)  •*2<-Y  1 ( I)  **2 
bO  AAA ( I , I )=0 . 5*( U/V) •SO( I) 

QO  32  K=1,NC 
00  30  J=1,NT 
30  CA(J,<) =CC( J,K) 

32  CONTINUE 


♦ * • 


C the  eigenfunctions  are  OETERMINED  OY  ITERATION.  XM 
C GIVES  THE  3AXIHUH  DEVIATION  OF  THE  COMPONENTS  OF  THE  CURRENT 
C KTH  ITtRATE  FROM  THOSE  OF  THE  PREVIOUS  K-1  TH  ITERATE. 

C the  IT-.RATIOi4  PROCESS  TERMINATES  IF  EITHER  K = 20  OR 
C XM  L.T.E.  1,*10**-3.  THt  EIGENFUNCTIONS 
C CORRES'^ONO  INC  TO  THE  EIGENVALUE  PI  OF  THE 
C OF  the  TRA  ISPOSE  OF  AAA  ARE  OETERMINED  FIRST. 


DO  140  K1=1,NC 
DO  120  <=1,20 
XM  = a 

00  100  J=l,NT 
Cl ( J) =0 
OJ  30  1=1, NT 

90  CI(J»  =CKJ)+CA  (I,Kl)  ♦AAA(I,J) 
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CI(J)=CI(J)/PI 

AdS (CA ( J ,K1 » -Cl ( J) ) 

Xi=AMAXl(XM,XMl) 

100  CONTIMUl 

00  110  J=1,MT 
110  CA{J,K1 )=CI( J) 

IF  (l.E-Oi-XM)  120,130,130 
120  CONIINUt 
130  XA-KKDsXi 
NA(<1) =< 

IhO  continue 

C THt  EIGENFUNCTIONS  CORRESPONOING  TO  PI  OF  THE 
C matrix  AAA  ARE  OcTERMINEU  NEXT. 

00  230  <1=1, NC 
DO  210  <=1,20 
XH=0 

DO  190  1=1, NT 
Cl ( I ) =0  . 

00  180  J=1,NT 

180  CI(I)=CI(I)fAAA(I,J) *CC( J,<1) 

CI(I»=CI(I)/PI 

XM1  = ABS (CC ( I ,<1) -Cl ( I)  ) 

XM=AMAXKXM,  XMl) 

190  CONTINUE 

00  200  1=1, NT 
200  CC(I,<1)=CI(I) 

IFd.E- 08-XH)  213,22  0,220 
210  CONTINUE 
220  XME(<1)=XM 
NI  (<1)  = < 

230  CONTINUE 

C THE  EIGENFUNCTIONS  HAVE  BEEN  OcTERMINEJ.  NI(Kl)  GIVES  THE 
C NUMBER  OF  ITERATIONS  IN  THE  DETERMINATION  OF  THE  <1TH  EIGEN- 
C function.  XME(<1)  GIVES  THE  MAXIMUM  DIFFERENCE  BETWEEN  THE 
C COMPONENTS  OF  THE  LAST  TWO  ITERATES  OF  THE  KITH  EIGENFUNCTION. 

RETURN 

END 


subroutine  PITER(N,<,£M,A,8,U 
DIMENSION  A ( 60 , 80) , 3 (80)  , V(80»  ,8C ( 30) 


c This  subroutine  solves  the  linear  system  of  equations 
C (PI*I-A)X=B  JY  AN  iteration  PROCESS.  HtRE  PI  IS  AN  EIGEN- 

C VALUE  OF  THE  MATRIX  A,  FIRST,  REWRITE  THE  SYSTEM  OF  EQUA- 

C TIONS  IN  THE  -ORM  ( I- A/PI ) X=8/PI  AND  SET  AC=A/PI  AND  BD=8/PI. 

c The  matrix  ac  has  i for  an  eigenvalue,  we  suppose  removed 
c from  the  vector  80  any  eigenvectors  corresponding  to  the 

C eigenvalue  1.  ONE  OBTAINS  X = BO  AC*  BO*' AC*  ( AC*  30)  ♦ ETC. 

C PITER  CALLS  THE  FOLLOWING  SUBPROGRAM! 
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o n 


MATV  THIS  3U3R3JTINE  GIVES  THE  PROQJOT  OF  A VECTOR  BY 

A matrix  amj  oivioes  the  resulting  vector  by  pi. 

**••♦♦**•*•♦ 

PI=i. 14 15R26535398 
00  5 1=  1,N 
S V(I»  =B( I) 

K = 1 

15  call  MATV( N, a, V) 

DO  20  I =1,N 
20  D ( I) =3<  I) fVC I) 

EM=  3 

DO  30  1=1, N 
U = 4 JS (V  (I)  ) 

E = AHAX1  (EM,U) 

30  EM=£ 

IF  (1.E-09-EH)  ^0,40, TO 

hO  IF(20-K)  60,60,50 
50  K=K+1 

GO  TO  15 
60  L = 0 

GO  TO  90 
70  L=1 

90  return 
ENO 


SJD^OUTINE  PMAT2(N,M,A,U, V> 
dimension  A(2,2),U(83,2),V(80,2» 


c This  subroutine  constructs  the  colFfil eent  matrix  used  in 

C OETERMlNIi'JG  THE  LINEAR  CJM3INATI0N  OF  EIGENVECTORS 
C CORRESPONDING  TO  THE  EIGENVALUE  pI  IN  A VECTOR. 


00  30  1 =1, M 
00  20  J=1,M 
DO  10  K=I,N 

10  S = S^  J (K  ,1)  *V  (K  , J» 

20  A(I,J)=S 

30  continue 

RETURN 
END 


SJ3ROUTINE  MATJ(N,A,V) 
dimension  h(80 ,80) ,V  (80)  ,S(80) 


C this  SU3ROUTIN--  HJLTIPLIES  A VECTOR  BY  A MATRIX  AND  THEN 
C DIVIDES  the  resulting  VECTOR  BY  PI. 


PI  = 5. 14  15926535998 
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DO  20  1=1, N 
S ( I ) =0 
DO  10  J=1,N 

10  sm  =S(  I)  ♦Ad,  j)»y/{  J) 
20  C0NIINU£ 

00  30  1 =1, N 
30  V(I)=S(I)/PI 
<ET 
END 


MTVlLH  ( I APE^  , INPUT  , OUTPUT) 

CONiON/  JK3/5K  8 0)  , 52(30), X(80),X1(83), 

2 X2(30)  ,Y(d0),Yl(80)  ,Y2(80),X1P(80), 

3 Y1P(80  ),TI(2)  , 50(30) 

CONMON/  -3K4/AA(  3 0)  ,03  (80)  ,CC(80,2)  ,CI  (80)  ,CA  ( 80,2) 
COMNON/ 1K5/AAA ( 80 ,30 ) , IP( 80 ) 

CONMON/ JK8/N (2 ) ,NC,NT 

CONI ON/ 8K7/NI(2),XN£(2),nA(2),XAM(2) 

C0NN0N/3K8/X0( 2) , YO ( 2) , N ( 2) , A ( 51 , 2) , 3 ( 40, 2) 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


NTVELP  IS  THE  CONTROL  PROGRAM  FOR  DETERMINING 
THE  tangential  VELOCITIES  INDUCED  BY  A DISTRI- 
BUTION OF  VORTICES  ON  THE  CONTOURS. 

IN  this  program  systems  of  linear  equations 

ARE  SOLVED  DIRECTLY, 

NC  NUM3ER  OF  CONTOURS 

N(I)  NUMJER  OF  SUBDIVISIONS  OF  THE  ITH  CONTOUR 

NT  TOTAL  NUMBER  OF  SUBDIVISIONS  ALL  CONTOURS 

NT=3UM  N (I)  , 1 = 1, NC 

HE  HAVE  USED  FIXED  RATHER  THAN  VARIABLE  DIMENSIONS.  IF  THE 
NUMBER  OF  CONTOURS  OR  IF  THE  NUMBER  OF  SUBDIVISIONS  IS  IN- 
CREASED THEN  THE  DIMENSION  STATEMENTS  HILL  HAVE  TO  BE  CHANGED 
AND  THE  PROGRAM  MTVELP  AND  SUBROUTINES  WHICH  IT  CALLS  WILL 
HAVE  TO  BE  recompiled  TO  REFLECT  THE  NEW  VALUES  OF  NT  AND  NC 
XO(K),YO(K)  APPROXIMATE  "CENTER"  OF  THE  KTH  CONTOUR 
A(L,K) ,3(L  ,K)  FOURIER  COSINE  AND  SINE  COEFFICIENTS  FOR  THE 

radius  VECTOR  R FROM  XJ(K),YO(K)  TO  A POINT  ON 
The  KTH  CONTOUR. 

THE  CENTERS  ( X J ( K ) , YO ( K ) ) AND  THE  FOURIER 
COEFFICIENTS  A(L,K)  AND  B(L,K)  ARE  READ  FROM 
THE  files  ,TAPt5.  THESE  QUANTITIES  MUST  BE 
PREPARED  AND  STORED  THERE  IN  ADVANCE  AND  IN 
THE  PROPER  FORMAT, 

AA,BB  CONTRIBUTIONS  TO  THE  TANGENTIAL  VELOCITY 

RESULTING  FROM  THE  VELOCITY  COMPONENTS 
AT  infinity. 

CC  COMPONENTS  0^  THE  EIGENVECTORS.  THESE 

PROVIDE  THE  CIRCULATION. 

CA  EIGENVECTORS  OF  THE  ADJOINT  PROBLEM 

OR  FOR  THE  TRANSPOSED  MATRIX. 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

u 

c 

c 

c 

c 


X(I),V(I)  COO-XDINATti  JP  Tmt  ith  STATION  ON  THE  CONTOJRS 

Xl(I),yi(I)  AKE  THE  EKST  OtRIVATIVES  W.R.T, 

Hl  central  ancle  theta  and  SO  ON 

XI P ( I) , Y 1P( I ) ARE  THE  FIRST  DERIVATIVES  W.R.T. 

A^C  lencth  and  so  on 

S1,S2  D::NJTE  the  1ST,  2ND, ETC.  DERIVATIVES 

OF  THE  ARC  LENGTH  W.R.T.  THE  CENTRAL  ANGLE 

SO  CHANGE  OF  INTEGRATION  VARIABLE  FACTOR 

HTVELP  calls  THE  FOLLOWING  SUBPROGRAMS 

mvalp  This  subroutine  uses  the  quantities  in  common  block 
a TO  COMPUTE  the  quantities, primarily  in  3K3. 

MP/IAT  THIS  SUBROUTINE  USES  THE  VALUES  IN  COMMON  BLOCK  3 TO 
COMPUTE  TH:  matrix  ErjTRIES  FOR  DETERMINING  THE 
tangential  VELOCITIES  AND  THE  E I GEN VEC T ORS . 

The  matrix  entries  are  stored  in  common  BLOCK  5. 

LUJ  A routine  for  the  DIRECT  SOLUTION 

OF  SYSTEMS  OF  LINEAR  EQUATIONS  TAKEN 
FROM  THE  ARL  LINEAR  ALGEBRA  LIBRARY, 

LU3  INCLUDES  LEQS3. 

THE  input  UATA  TO  RUN  THIS  PROGRAM  IS  AS  FOLLOWS 

*EOR 

B3BNC 

B 19N ( 1 ) 

ETC. 

JBBN(NC) 

*£OR 


NAMElIS  T/NAM/l  LL/NAM  1/LLL/NAM2/LLL/NAM3/LLL/NAM*+/LLL/NAM5/ 
2 LLL/NAH6/LLL/ N A M7/ L L L / N A M 3/ L L L/ NAM9/ L LL 
LLL  = 1 

PI= i. 14 15B26535aBtt 
READ  90J,NC 
NT  = 3 

00  10  I=1,NC 
READ  90  0, N(I) 

10  NT=NT*N  (I) 

READ (5 , NAM4) 

DO  20  J=1,NC 
REAJ(5,  UOQ)M(J) 

MP=M( J) +1 
MM=  M ( U)  -1 

READ (b,  901 ) XO( U)  , YD ( U) 

READ(5,901)(A(I,U),I=1,MP) 

20  REA  1 ( 5,  901)  (8( I , U) , 1 = 1, MM) 

CALL  MVALP 
call  MPiMAT 

call  lu 3(nt ,aa a ,ao, IP) 
call  LE QS3 (NT, a AA, 50 ,IP, AA) 
call  LEQS3 (NT, AAA, 30 , IP, BB) 

R^AD (b, NAM5) 

06  -*0  Kl=l,NC 
11=1 
12=0 

WRITE (5 ,9)0) NA ( Kl> 
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w^ari  (&  ,90  1)  xAi  (Ki) 

00  5 0 K=1,NC 
I^  = I^♦N  (K) 

W(?irE(5,931)  <:A(I,<1),I  = I1,I2) 
M-?IT£  (5  ,900)  N(  <) 

30  I1=I1+N(K) 

4C  CONTINUE 

Wf^ITE  (5  ,NAM6) 

11  = 1 
12  = 0 

00  50  <=1,NC 
I2=I2  + N (K) 

WRITE (5 ,90 1) (AA(I),I=I1,I2) 
WRITE<5  ,900)N(K) 

50  I1=I1+N(K) 

11  = 1 
12=0 

00  50  K=1,NC 
I2=I2  + N <K) 

WRITE(5,901) {ja(I),I=Il,I2) 
WRITE (5  ,900) N( K) 

60  I1=I1*N(K) 

DO  tiO  K1=1,NC 
11  = 1 
12=0 

WRITE (5  ,90  0) NI  ( Kl) 

WRIT£(5,901)  XNECKl) 

03  70  K=1,NC 
I2  = I2tN  (K) 

WRITE(5,901) (CC(I,<1) ,1=11,12) 
WRITE(5,900) N(K) 

70  I1=I1+N(K) 

80  CONTINUE 

STOP 

900  F0R1AT(I5) 

901  FOR-IAT  ( 5E23.  1h) 

ENO 


SU9ROUTINE  NPNAT 

CO-ITON/ 3K3/j1(  80)  , 32(90), X(30),X1(80), 

2 X2(90)  ,Y(90),''1(80),Y2(30),X13(80), 

3 Y1P(90  ),TI(2)  ,SO(90) 

C0110N/ 9K4/AA(rtO) , 33 ( 8 0 ) , CC ( 8 0 , 2 ) , C I ( 8 0) , C A { 80 , 2) 
COTTON/  3«;5/AAA  (d0,30),IP(60) 

C01NON/ JKb/N(2) ,NC,NT 

CJHION/  )K7/NI(2),XNE(2),NA(2),XAM(2) 


C THIS  SU3R0UTINE  DETERHINES  TmE  MATRIX  ILEMENTS  NEEDED  FOR  THE 
C COMPUTATION  OF  THE  TANGENTIAL  VELOCITIES.  THE  MATRIX  FOR 
C COMPUTING  THE  EIGENFUNCTIONS  IS  DETERMINED  FIRST. 
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«*»*»**»»»»«» 


PI  = i. m 15^265353  )6 

i2=Nr-i 

TO  ?0  1=1,12 

OQ  20  J=J1,NT 

U=(<(I)-X(J))*V1P(I)-(Y(I)-Y(J))*X1P(I) 

V=(X(I>-X( J))**2*(Y(I)-Y(J))**2 

20  a^a ( I , j ) = (u/ V) *so ( J) 

30  CONTINUE 

DO  ^0  1=2, NT 
J2=I-1 

00  +0  J=l,  12 

U=(<(1)-X(J))*Y1P(I)-(Y(I)-Y(J))*X1P(I) 

V=(X(I)-X(J»)**2+(Y(I»-Y(J))**2 

40  au(i,j)  = (u/v)  »sj(j) 

50  CONTINUE 

DO  -jO  I =1,  NT 

■ U=Y2( I) ’XlPi I) -X2(I) *Y1P(I) 

V=X1 (i ) »*2^Y 1 ( I )**? 

60  aaa(i,i)=o.5*(j/v)*so(i) 

♦ *****»»*******»*******»**»)f****»***4***4*****4*******»****»»-* 

C THE  lIjENFUNCTIONS  ARE  OETtRMINED  OY  ITERATION.  XM 
C GIOES  THi  HaxiHUM  DEVIATION  OF  THE  COMPONENTS  OF  THE  CURRENT 
C KTH  IT-.RAT.  FROM  THOSE  OF  THt  PREVIOUS  K-1  TH  ITERATE. 

C THE  IT-RATION  PROCESS  TERMINATES  IF  EITHER  K=20  OR 
C XM  L.T.E.l  .*10**-3.  AFTER  THE  EIGENFUNCTIONS  ARE  JETERMInEO 
C THt.  MATRIX  MOST  3E  MOOIFIEJ  TO  DETERMINE  THE  CONTRIBUTIONS 
C TO  THE  tangential  VELOCITY  DUE  TO  A UNIT  X ANU  Y COMPONENT 
C OF  THE  velocity  AT  INFINITY. 

C COMPUT.  THE  EIGENVECTORS  CO TRESPONOING  TO  THE 
C eigenvalue  PI  FOR  THE  AOJOINT  PROBLEM. 


DO  62  K=1,NC 

00  60  J=1,NT 

60 

CA( J,K) =CC(J,<) 

X a M ( < ) = 0 

Nt (<) =0 

82 

CONTINUE 

C COMPUTE  THE  EIGENVECTORS  CORRESPONDING  TO  THE 
C lIGENVALUE  pi  for  the  matrix  AAA. 


00  230  <1=1, NO 
00  210  <=1,20 
XM=0 

□0  IRO  1=1, NT 
Cl ( I ) =0  . 

00  130  J=1,NT 

160  CI(I)=CI<IM-aAA(I,J»  *CC(J,K1) 
Cl  ( I ) =C  Id  ) /PI 
XM1=A6S (CC (I ,K1 » -Cl ( I) ) 
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XM  = 4iAX  KXM,  XMl  ) 

190  CONTIMU£ 

03  200  1=1, ST 
200  CC( I ,K1  )=CI ( I) 

IF(l.L-08-XS)  210,220,220 

210  continue: 

220  XS£(K1)=XM 
NI (Kl) = < 

230  CONTINUt 


C THt  EIGENF.JNCTIOSS  HA\/£  BEEN  DETERMINED.  NI(Kl)  GIVES  THE 
C NUMBER  OF  ITERATIONS  IN  THE  OETERMINATI ON  OF  THE  KITH  EIGEN- 
C FUNCTION.  XME(Kl)  GIVES  THE  MAXIMUM  DIFFERENCE  BETWEEN  THE 
C COMPONtSTS  OF  THi  LAST  TWO  ITERATES  OF  THE  KITH  EIGENFUNCTION. 


c the  matrix  elements  for  the  modified 

C MATRIX  ARE  COMPUTED  NEXT. 


DO  310  1=1, NT 
03  JOO  J=1,ST 
U=0 

DO  290  <=1,NC 
290  U = U«-CC(  I,K»*CA  ( J,K) 

30  0 AAA  ( I , J )=-AAA  ( I , JH-PI*U*SD(  J) 

310  CONTINUE 

03  320  1=1, NT 

320  AAA(I,I  )=PI*-AAA  (1,1) 

RETURN 

900  FOR-1AT  ( •♦15) 

901  FORMAT  ( 5E23.  1**) 

END 


PKO  MT\/£LS(TaP£5, inpot,  OUTPUT) 

COHN  ON/  -3T1  /aa  ( d J ) , ) 3 (80  ) ,CC  (3  0,2) 

2 /3T2/TI(2),SJ(3J),S1(30)  ,32(30)  , 33(30), X(6  0 ), 

3 Xl(30)  ,X2(  3 0)  ,X3(30),Y(30),Yl(30),y2(80),Y3(30),XlP{dO), 

4 x2°(30),X3P(8d)  ,Y1P(30)  ,Y2P(80),Y3P(80) 

5 /3T3/DAa(80) ,03B(30) ,OCC(80,2) 

D /3T4/C(2),XX(2),YY(2),M{2),A(51,2),B(49,2) 

7 /9T5/A AA ( 33 , « 0 ) , IP ( 80 ) 

3 /3T5/TaA(30),T)0(80) ,TCC(80,2) 

9 /3T7/N  (?)  ,NC, NT 

NAN£LIST/NA1/LLL/NANl/LLL/NA'12/LLL/NAM3/LLL/NAN4/LLL/NaN5/ 
2 LLL/NA  N6/LLL/  NA,N7/LLL/NAM3/LLL/NAH9/LLL 


NC 

N (I) 
NT 


W£  HAY; 


NUN‘3Er<  OF  CONTOURS 

NJN'JtR  OF  SUBDIVISIONS  OF  THE  ITH  CONTOUR 
total  number  of  SUBDIVISIONS  ALL  CONTOURS 
NT  = SUM  N (I)  , 1 = 1, NC 

USED  FIXlD  rather  THAN  VARIABLE  DIMENSIONS.  IF  THE 
NUMBER  OF  CONTOURS  OR  IF  THE  NUMBER  OF  SUBDIVISIONS  IS  IN- 
CREASED th.n  The  dimension  siatements  /Ull  have  to  be  changed 
AND  The  program  MTVELS  AND  SUBROUTINES  WHICH  IT  CALLS  WILL 
HAVE  TO  BE  recompiled  TO  REFLECT  THE  NEW  VALUES  OF  NT  AND  NC 
XX(K),YY(K()  approximate  "CENTER”  OF  THE  KTH  CONTOUR 
A (L  ,<)  , 3 (L  ,K)  FOURIER  COSINE  AND  SINE  COEFFICIENTS  FOR  THE 

RADIUS  VECTOR  R FROM  XX(K),YY(K)  TO  A POINT  ON 
The  KTH  CONTOUR. 

The  cenTcRS  ( xx (k>  , yy (k ) ) and  thl  fourier 
COEFFICIENrS  A(L,<)  AND  3(L,K)  ARE  READ 
FROM  The  file,  tapes,  these  ouantities 

MUST  Bl  prepared  AND  STORED  IN  ADVANCE 
IN  THE  PROPER  FORMAT. 

AA,83,C:  PICTITIOUS  DENSITY  FUNCTIONS  RESULTING  FROM  THE 

velocity  components  at  infinity  and  CIRCULATION 
ABOUT  EACH  CONTOUR 

OlRIVATIVES  of  the  DENSITY  FUNCTIONS  W.R.T.  THE 

integration  variable  theta 
the  contributions  to  the  tangential  velocity 

COORDINATES  OF  THE  ITH  STATION  ON  THE  CONTOURS 

xi(i),ri(i)  ARE  the  first  Derivatives  w.r.t.  th 
C-NTRAL  angle  theta  and  so  on 
XIP(I) ,Y1P(I)  are  the  first  derivatives  w.r.t. 
ARC  LENGTH  AND  SO  ON 

denote  the  1ST, 2ND  E TC , DE R I VA T I V ES  OF  THE 
ARC  length  W.r.t.  the  central  angle 
change  OF  INTEGRATION  VARIABLE  FACTOR 
MTVELS  CALLS  THE  FOLLOWING  SUBPROGRAMS 

MVALS  IHIS  SUBROUTINE  USES  THE  QUANTITIES  IN  COMMON  BLOCK  4 
TO  COMPUfiL  THE  QUANTITIES  IN  COMMON  BLOCKS  1 AND  2 


OAA , DDB, OCC 

TAA , T3  J, TCC 
X (I ) , Y ( I ) 


S 1 , S 2 , S 3 


SO 


90 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


MOENS  THIS  SUB-<OUT  US^S  THE  VALUES  IN  COMMON  BLOCK  2 TO 
COMPUTE  THE  MATHIK  ENTPILS  FUP  DETERMINING  THE  OEN-. 
SITIES.  THE  MATRIX  ENTRIES  ARE  STORED  COMMON  BLOCK  5. 

LU3  A ROUTIN’  FOR  THE  DIRECT  SOLUTION 

OF  S/STEMS  OF  LINEAR  EQUATIONS  TAKcN 
FROM  The  arl  linear  algebra  library, 

LEQS3  IS  II4CLUDED  IN  LU3 

MODEMS  THIS  SUBROUTINE  COMPUTES  THE  DERIVATIVES  OF  THE  DENSITIES 
MTVS  this  subroutine  COMPUTES  THE  CONTRIBUTIONS  TAA,T3B  ANO  TCC 

TO  the  TANGENTIAL  VELOCITIES. 

THE  I4RUT  DATA  TO  RUN  THIS  PROGRAM  IS  AS  FOLLOWS 
NOTE-THE  LETTER  B DENOTES  A BLANK  IN  THE  INPUT  DATA. 

*£OR 

BBBNC 

BBBN(l) 

ETC. 

33BN(NC) 

♦EOR 


LLL=  1 

PI=i. 14159265353970 
READ  90  0, NC 
NT  = Q 

DO  10  1=1, NC 
READ  900, N(I) 

10  NT=NT+N(I) 

REAJ(5, NAM4) 

DO  20  U=1,NC 
READ(5, 900)M(J> 

M1P=M  ( J ) 4-1 
M1M=M ( J )-l 

R£AD(5,901)XX(U),YY(J) 

READ (5, 901)(A(I,J),I=1,M1P) 

20  read <5,  BOl  ) (B( I ,U» , I =1, MIM) 

call  mvals 
CALL  NDENS 

CALL  LU3(NT,AAA,30,IP) 

CALL  LE  IS3 (NT , AAA,fl0 , IP, AA) 

CALL  LEQS3(.NT,  AAA,90  ,IP,  3B) 

DO  25  J=1,NC 

25  call  LE 3S3 (NT, AAA ,60 , IP,CC  (1, J) > 
R£A0(5, NAM5) 

WRir£(5,900)  NT 
WRITE(5,901) (AA(I),I=1,NT) 

WRITE (5  ,90  0) NT 
WRITE(5,901( (BJ(I),I=1,NT) 

DO  3J  J=1,NC 
WRITE (5  ,900) N( J) 

30  WRITE (5  ,90 1)  (CC ( I,U)  ,I  = 1,NT) 
WRITE(5  ,NAM6) 

call  modems 

WRITE (5 ,900) NT 
WRITE(5,90l)  (0AA(I),I  = 1,NT) 


An-A042  984 


UNCLASSIFIED 

2of  2 
^%42S»l 


AIR  FORCE  FLIGHT  DYNAMICS  LAB  WRIGhT-PATTERSON  AFB  OHIO  F/B  20/4 
AW^77*^C^L^KEL^R^^^^*^°^  dimensional  INCOMPRESSIBLE  FL— ETC(U) 

AFF0L-TP-77-a7  « 


WRITE (5 ,900) NT 
WRITE{5,901) (039 (I), 1=1, NT) 

00  ^0  J=1,NC 
WRITE (5 ,900) N( J) 

40  WRITE(5  ,901)  (OCC(I,J)  ,I  = 1,NT) 
WRITE  <5  ,NA97 ) 
call  MT  \/S 
WRITE (5 ,900) NT 
WRITE(5,901)  (TAA(I),I  = 1,NT) 
WRITE(5  ,900)  NT 
WRITE(5,901)  (T9  3(1), 1 = 1, NT) 

00  50  J=1,NC 
WRITE (5  ,90  0) N( J) 

50  WRITE(5  ,901)  (TCCd,  J)  ,I  = 1,NT) 
WRITE (5 ,NA99 ) 

STOP 

900  F0R9AT(I5) 

901  FORMAT ( 5E23. 14) 

ENO 


SU9R0UT  INE  MVALS 

COMMON/  TTl/AA(aO>,a3(dO),CC(dO,2) 

2 /9r2/TI(2),SO(30),Sl(80),S2(80»,S3(90),X(d0), 

3 Xl(30)  ,X2  (90),X3(8a),Y(80),Vl(80),Y2(80),Y3(80),XlP(80), 
H X2P(80),X3P(80) , YIP (80) ,Y2P(80» ,Y3P(80) 

6 /3T4/C  (2)  ,XX(2)  ,YY(2)  ,M(2)  ,A(51,2),0(49,2) 

9 /BT//N  (2)  ,NC, NT 
REAL  MFOU,MJF0U,MO2FOU,MO3FOU 


C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


this  SJ9R01JTIN_  prepares  tables  of  VALJES  TO  BE  USED 
IN  LATER  COMPUTATIONS.  SPECIFICALLY,  IT  COMPUTES  THE 
C00R0INATE3  AT  STATIONS  ALONG  THE  CONTOUR  C,  THE 
OERIVATIVES  OF  THESE  COORDINATES  W.R.T.  THE 
POLAR  7ARIA3uE  THETA=T,  AND  WITH  RESPECT  TO  ARC 
length,  these  are  DENOTED  E.G.,  9Y  XI  AND  XlP  RESP, 
IT  ALSO  COMPJT-S  THE  DERIVATIVES  OF  THE  ARC  LENGTH 
W.R.T.  THE  POLAR  variable.  THESE  DERIVATIVES  ARE 
UENOTEJ  UY  S1,S2,  ETC. 

THIS  SUBROUTINE  USES  THE  FOLLOWING  RcAL 
FUNCTION  PROGRAMS. 

MFOU  lVAlJATES  a trigonometric  POLYNOMIAL  WITH 
FOURIER  coefficients  A(L,K)  AND  B(L,K)  AT 
A GIVEN  ANGLE  (IN  RADIANS). 

MDFOU  THE  DlRIVATIVE  OF  MFOU  W.R.T.  THETA. 

M02F0U  THE  DlRIVATIVE  OF  HOFOU  W.R.T.  THETA, 

M03FOU  THE  DERIVATIVE  OF  MD2FOU  W.R.T,  THETA. 


PI=3. 14 159265358898 
T = 0. 

11  = 1 
IL  = 0 
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DO  35  J=1,N0 
IL=IL>N  (J) 

TKJ)  =2.*PI/N(  J) 

00  30  I=II,IL 
«=HP0U( T,^ (J) , ft, 3, J> 

R1=1DF0U(T (J) ,ft,3, J) 

R2=iJ2F0U(T, M( J) ,ft,9,J) 

R3  = '1J3F0U(T,M(J),ft,0,J) 

SI  ( I)  =S  JRT  (R**24'R1**2) 

SD(I»  =S1(I»*TI  ( J) 

U=?+R2 

V=R1*U 

S2(I>  =V/S1 (I ) 

W = R2*L)^R1*  (RltR3) 

02=S1(I)**2 

S3 ( I) = ( 31( I) *W-V*S2 ( I) » /02 
C0=C0S(T) 

SI=3IN( T) 

X(I)=XX (J) +CO*R 
Y(I»=YY (J) 

XKI)  =C0*R1-R»SI 

Y1(I)=SI*R14-R*00 

X2( I>  =C0*R2-2*SI*R1-C0*R 

Y2(I)  =S  I*R2«-2*CO*Rl-SI*R 

X3(I) =R3*C0-3. *R2*SI-3.*R1*C0^SI*R 

Y3(I) =R3*SI^3. *R2*C0-3.*R1*SI-C0*R 

03=31 (I )*D2 

05=02*03 

X1P(I)=X1{I)/S1 (I) 

Y1P(I)=Y1(I)/S1(I) 

U=S1(I) *X2(I)-X1(I)*S2(I) 

X2P ( I) = J/03 

V = S1(I) *r2 (I )-Y 1(1) *S2  (I) 

Y2P(I)  = \//D3 

Z=SI<I)*(31(I»*X3(I)-X1(I)*S3(I)) 

W = S1(I»*(S1(I)*Y3(I»-Y1(I»*S3(D) 

X3P(I)= (Z-3.*U*S2(I) )/05 
Y3P(I)=  (M-3.*V*S2(in/D5 
Aft(I)=YlP(I) 

BB<I)  =X1P(I) 

00  20  K=1,NC 

U=(XX{K»-X(I»>*X1P(I)*(YY(K)-Y(I))*Y1P(I» 

U=-iJ 

V=(XX(K)-X(I))**2*(YY(K)-Y(I))**2 
20  CC(I,K)=U/7 

T = T + TI(  J) 

30  CONTINUE 

II=II*N  (J) 

T = 0. 

35  continue 

RETURN 
ENO 
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SJ30UTIN£  MOwNS 

CO'HON/  JT2/TI  ( 2>  ,SJ(3G),S1(90),S2(J3),S3(S0),X(30), 

3 Xl(30),X2(30»,xndO),Y(8C).Vl(80),Y2(80),Y3{dO),XlP<80) 

4 X2»(80),X3P(90),Y1P(9G>,Y2P(90),Y3P(80) 

6 /3T4/C (2),XX(2),YY(2),M(2),A<51,2),B(49,2) 

7 /3T5/A 4A( 90, 8G ) ,IP(80) 

9 /9T7/N  (2)  ,NC, NT 

C THIS  SUBROUTINE  30NPUTSS  THE  COEFFICIENT  MATRIX  FOR 
C DETERMINING  THE  SOURCE  DENSITIES. 


PI=3. 14 159265359979 
I2=NT-1 
DO  30  1=1,12 
J1=I*1 

00  20  J=J1,NT 

U=(X(I)-X(J))*Y1P(I>-(Y(I)-Y(J))»XIP(I) 
Y=(X{I)-X(J))**2+(Y(I»-Y(J)»**2 
20  AAA (I , J )= ( U7 Y) "SOC J) 

30  CONTINUE 

DO  50  1=2, NT 
J2=I-1 

00  ‘♦O  J=1,J2 

U=(X(I»-X(J))*Y10(II-(Y(I)-y<J))*XlP(I) 
V=(X(I)-X(J))**2«-(Y(I)-Y(J))**2 
40  AAA(I, J)=(U7V) *30(0) 

50  CONTINUE 

DO  60  1=1, NT 

U=Y2(l)*XlPlI)-X2lI>»YlPm 
V=X1 ( I)  *»2+Yl( I>  **2 
60  AAA(I,I )=PI+0. 5*(U/y)*SO(I) 

RETURN 

ENO 


SUBROUTINE  MDOENS 
DIMENSION  SM3(2) 

CO 1MON/ 3T1/AA( 93) , 93 ( 80 ) , CC ( 3 0 , 2) 

2 79T2/TI(2),SJ  (30), SK30), 32(80), S3(80),X<80), 

3 XK80)  ,X2  (30),X3(90),Y(80),Yl(80),Y2(80),Y3(80),XiP(80) 

4 X20(90),X3P(90),Y1P(30),Y2P(80),Y3P(80) 

5 /9T37JAA(90),J33(80),OCC(80,2) 

6 /3T47C(2),XX(2),YY(2),M(2),A(51,2),B(49,2) 

9 79T7/N  (2) ,NC, NT 


C THIS  SUBROUTINE  COMPUTES  THE  OERIVATVES  OF  THE 
C JENSITIY  FUNCTIONS  W.R.T.  ThE  INTEGRATION 

C VARIABLE  T (=  THETA) 

PI=3. 14 15926535999 

SM1=3M2=0 

00  5 K=1,NC 
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5 S‘13{<)  = 0 

1 = 1 

GD  TO  30 
10  SMi=SM2=0 

00  15  K=1,NC 
15  313(0  = 0 

J2=I-1 

00  20  J=1,J2 
U=X(I) -X(J) 

>/  = Y{  I)  - Y(J) 

Z=U» JVV*V 

W = z* (U*Y2? (I )-V*X2P{  I) ) 

W = K-2.*  (U*Y1P(I)-V*X1P(I)  )*(U*X1P(IM-V*V1P(I)) 
W= (1*S0  ( J) ) / (Z»Z) 

SMl=SMl+W*Aft (J) 

312=312 +M*33(J) 

00  17  K=1,NC 

17  313(0  = SM3(K)«-W*GC(  J»K) 

20  CONTINUt 

IF(I-NT)  30,50,50 
30  Jl  = I«’l 

00  40  J=J1,NT 
0=X(I)-X(J) 

Y=Y(I)-Y(J) 

z=u*u*v*v 

W = Z*(0*Y2P(I)-V*X2P(D) 

M = 1-2.  ♦ (U*  Y1P(  H -V*X1P(  I)  ) ♦ (U*X1P(  I)  <-V*YlP(  I)  > 
W=(1*S0  (J) »/ (Z*Z) 

311=311VM*AA(J> 

312=312  03  (J) 

00  55  K=1,N3 

35  S13(0=SM3(0  + W*CC(J,0 

40  CONTINUE 

50  Z=X1P( I ) **2»Y1P ( I) **2 

W1  = X1P( I»*Y  3P( I )-YlP (I>*X3P (I ) 

U=X1P (I »*Y2P (I J -Y1P( I) *X2P( I) 
y/  = XlP(I»*X2P(IM-YlP(I)*Y2P(I) 

W = S3(I)  ♦(2.*Z*W1-3.*U*«0/(6.*Z*Z) 


C THE  VALUES  OP  THE  INTEGRALS  ARE 


Sli=(SMl^W*AA(I)) 

312= {S12»H*BB(  I)  I 
00  52  K=1,N0 

52  313(0  = (S13(K)  «-W*C:  (1,0  ) 

0=31(1) **2 

Hl=(31(I)*Y2(I)-Yl(I)*S2(I))/0 
W2=(S1( I)*X2(I) -X1(I)*S2(I) )/0 
OAA(I)= (W1-S1(I)*S11)/PI 
033(1)= (M2-S1(I)»S12)/PI 
00  55  K=1,NC 

Z=(X(I) -XX (O ) •*2^(Y (I) -YY (K) ) **2 

Z1  = 2.*((X(I)-XX(0)*X1(I)^(Y{I)-YY(0)*Y1(I)> 

U=  (Y  (I)  -YY  (O  ) • Yl’3(  I)  ♦(X  ( I)  -XX(K) ) ♦XIP(I) 
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V=(y(I)-Yr(<))*(Sl(I)*Y2(I)-Yl(I)*S2(I)) 

V=\/*(X(I»-X)<<K)»*(S1(I)*X2(I)-X1(I)*S2(I)) 

U1=1»//D)*Y1(I)*Y1P(I)+-X1(I)*X1P(I) 

Wl= ( / (Z*Z) 

53  DCC(I,K)=(W1-S1(I)*SM3(K) )/PI 

1  = 1+1 

IF  (I-NT)  10,10,60 
60  RETURN 
ENO 


subroutine  MTYS 
DIRENSIUN  SR3(2) ,W3(2> 

COiMON/3Tl/AA(ttO>  ,B3(80),CC(80,2) 

2 /3T2/TI(2),Sl)(80),S1(80>,S2(80),S3(80),X(80), 

3 XK80)  ,X2  (80),X3(30),Y(30>,Y1{80),Y2(80),Y3(80),X1P'80», 

4 X2P(80),X3P(30»,Y1P(80),Y2P(80),Y3?(80) 

5 /aT3/Oafl(3  0)  ,D.3B(30)  , 000(80, 2) 

b /3T-4/'0(?),XX(2),YY(2),M(2),a(51,2),3(49,2> 

3 /3T6/T AA( 30 ) , T 33(30  ) ,TCO (30 , 2) 

9 /3TZ/N  (2) ,NO,NT 
PI  = 3.  14  1532o5353979 

0 THIS  subroutine  OOHPUTES  the  CONTRIBUTIONS  TAA,T3B 
0 ANO  TOO  TO  THE  TANGENTIAL  VELOCITIES  AT  STATIONS 
0 AROUN3  THE  CONTOURS, 

C FIRST  COMPUTE  THE  VALUE  OF  THE  INTEGRANO  FI  AT  THE 

C STATIONS  U NOT  EOUAL  TO  I,  THE  COMPUTATION  OF 
C FI  SEEMS  TO  SEPARATE  NATURALLY  INTO  3 CASESI  1=1 
C 1 L.T.  I L.T.  N ANJ  I=N.  the  COMPUTATION  OF  FI 
C IS  FOLLOHEJ  3Y  the  evaluation  of  THE  INTEGRAL. 

1 = 1 

SM1  = SM2  =0 
00  13  K=1,N3 
10  SM3(<)=3 

00  50  U=2,NT 

U=  (X  ( I)  -X(  U)  ) ♦ X1P(I»  + (Y  (I)-Y(un  *Y1P  (I) 

V=(X(I)-X(J))**2+(Y(I)-Y(U))**2 

FI=iO(U)*(U/V) 

SM1  = SM1  <-AA  (U>  *FI 
SM2=SM2+33(U>*FI 
DO  20  K=1,NC 

23  SM3 (<) =SM3 (<> +CC (U,<» *FI 

50  CONTINUE 
GO  TO  250 
60  SM1=SM2=0 

00  o2  X=1,N0 
62  SM3(<>=3 

J2=I-1 

00  ^0  U=1,U2 

U= (X (I ) -X( U) )*X 1P(I) + (Y (I )-Y(U) ) *Y1P (I ) 
V=(X(U-X(U))**2+(Y(I>-Y(U)>**2 
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FI=30(J)*(U/V) 

S11=SH1 ♦AA ( J) *FI 
SH2=  312  ♦-B-3  ( J)  *FI 
00  o9  K=1,N3 

69  S13 (K) =S13 (<) ♦OC ( J,K) *FI 

70  continue 
J1=I+1 

00  90  J=J1,1T 

U=(X(I)-X(J»)*XlP(I»f(Y(I)-Y(J>)*Yl>3(I> 
V=(X(I)-X(J))*»2MY(I)-Y(J))**2 
FI=30(J)*( J/V) 

SH1=SM1 *AA (J)*FI 
S12=S12>B3 ( J)*FI 
DO  79  K=1,NC 

78  SM3(K)  = 3M3(<)4-0G(J,<)*FI 

90  CONTINUi 
GO  TO  250 
90  S11=SM2=0 

DO  92  K=1,NC 
92  S13(K)=0 

J2  = m-1 
DO  100  J=1,J2 

U=(X(I)-X(J))*X1P(I)+(Y(I)-Y(J))*Y1P(I) 
V=<X(I»-X(J))**2+(Y(I)-Y(J))»*2 
FI=SO( J)*(U/V) 

SM1=S11+AA ( J»*FI 
SH2=S12+B9(J)*FI 
00  99  K=1,N0 

99  S13(K)>=iM3(K)^CC(J,K)*FI 

100  CONTINUE 
GO  TO  250 

c C01PJTE  Integrand  at  ith  station  here. 

250  Z=X1(I) •*2*Y1(I)*»2 

U=(X1(I)*X1P(I)*Y1(I)*Y1P(I»)/Z 

V=Xl(I»*X2(I)fVl(I)*Y2(I> 

H = 0.5*(X1P(I)*X2(IM-Y1P(I)*Y2(I)) 

W1  = -U*{S1(I)  ’BAAdM-AAd)  ♦SZd)  ) 

W2=-U* { 31(1) *093(1) ♦8B(I) *S2 (I) > 
W1=W1*((AA(I)*31(I)) *(U*V-W)/Z) 

W2  = H2*(  (03(I)*S1(I)  ) ’(U^V-vO/Z) 

DO  260  <=1,NC 

M3  0O  =- J*(S1  (I)  *OCC(  1,0  ♦■CCdfK)  *S2(  I)  ) 
260  M3(K)=W3(OM(CC(I,<)*3l(I))*(U*V-W)/Z) 

C the  YALUES  of  THi  INTEGRALS  ARE  GIVEN  3Y 
Srtl=SMl>Wl*SO( I) /SI (I) 

S12  = SM24-W2*S0(I)/S1(I) 

DO  270  <=1,NC 

270  S13(K)  = 3H3(0*H3(K)*SO(I)/S1(I) 

C THE  TANGENTIAL  VELOCITIES  ARE  THEN 
DO  290  <=1,NC 

U=(Xd)  -XX  (O)  •Y1P(I)-(Y(I)-YY  (O  )*X1P(I) 
V=  (X  (I ) -XX  (O  ) ♦•2*(Y  (I  )-YV  (K)  ) **2 
290  TCC(I,K)=(U/V) *SH3(K) 

TAAd)  = X1P(I)-3M1 


T3H(I)=YlP(I)tSM2 
I = I«-1 

IF(NT-I)  J00,'JO,60 
300  RITURN 
901  FO'^'iAT  ( 5E2  3. 1h) 

END 


P-tJjRAM  NTVLSI  (TAPE5  , INPUT, OUTPUT) 

COMMON/  )T1/AA( 30) , 39 ( d 0 ) , CC ( 8 0 , 2) 

2 /3T2/TI(2),SJ(30),S1(30),S2(80),S3(90),X<80), 

3 X1(3C)  ,X2(30),X3(a3),Y(8C),Yl(80),Y2(80),Y3(90),XlP(80), 

4 X2P (80), X3P (90), Y1P(30)  ,Y2P(80)  , Y3P  (80) 

5 /9T3/DAA(90) ,D38(80) ,OCC(80,2) 

b /3T4/C(2),XX(2),YY(2),M(2),A(51,2),B(49,2) 

7 /3T5/AAA(90,90) ,I»(80) 
a /ITb/T  AA ( 80 ) , T60(8O ) ,TCC  (80, 2) 

9 /9T7/N  (2)  ,NC, NT 

NANELIS  T/NANYLLL/NAN l/LLL/NAi2/LLL/NAM3/LLL/NAN4/LLL/NAH5/ 
2 LLL/NA  16YLLL/NA'17/LLL/NAM9/LLL/NA89/LLL 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


NTVLSI  IS  the  CONThiOL  P'^OGRAM  FOR  DETERMIN- 
ING THE  TANGENTIAL  VELOCITIES  DUE  TO  A 
JISTRIJJTION  OF  SOURCES  ON  THE  CONTOURS. 

IN  THIS  PROGRAM  THE  LINEAR  SYSTEMS  OF  EQUATONS 
ARE  solved  by  an  ITERATIVE  PROCESS. 

NC  NUMBER  OF  CONTOURS 

N(I)  NUMBER  OF  SUBDIVISIONS  OF  THE  ITH  CONTOUR 

NT  total  number  OF  SUBDIVISIONS  ALL  CONTOURS 

NT  = SUM  N ( I)  , 1=1, NC 

WE  HAVE  USED  FIXED  RATHER  THAN  VARIABLE  DIMENSIONS.  IF  ThZ 
NUMBER  OF  CONTOURS  OR  IF  THE  NUMBER  OF  SUBDIVISIONS  IS  IN- 
CREASED THiN  the  dimension  statements  will  have  to  be  changed 
AND  The  program  mtvels  and  subroutines  which  it  calls  will 
HAVE  to  be  recompiled  TO  REFLECT  THE  NEW  VALUES  OF  NT  AND  NC 
XX(K),YY(K)  APPROXIMATE  "CENTER"  OF  THE  KTH  CONTOUR 
A (L  ,K)  , 3(L  ,K)  FOURIER  COSINE  AND  SINE  COEFFICIENTS  FOR  THE 

RADIUS  VECTOR  R FROM  XX(K),YY(K)  TO  A POINT  ON 
The  KTH  CONTOUR. 


The  centers  ( x<  (o  ,yy(o  ) and  the  Fourier 

COEFFICIENTS  A(L,K)  AND  B(L,K)  ARE  READ  FROM 

THE  file,  tapes.  these  QUANTITIES  MUST  BE  PRtPAREO 

AND  STORED  IN  ADVANCE  IN  THE  PROPER  FORMAT. 


AA, B3,CC 


DAA,033,DCC 

TAA,T33,TCC 
X (I  ) , Y ( I ) 


FICTITIOUS  DENSITY  FUNCTIONS  RESULTING  FROM  THE 
VELOCITY  COMPONENTS  AT  INFINITY  AND  CIRCULATION 
ABOUT  EACH  CONTOUR 

DERIVATIVES  OF  THE  DENSITY  FUNCTIONS  W.R.T.  THE 
INTEGRATION  VARIABLE  THETA 

The  CONTRIBUTIONS  TO  THE  TANGENTIAL  VELOCITY 
COORDINATES  OF  THE  ITH  STATION  ON  THE  CONTOURS 
X1(I),Y1(I)  ARE  THE  FIRST  DERIVATIVES  W.R.T.  TH 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


central  angle  theta  anj  so  on 

XIP(I)  ,Y1P(I)  APE  THE  FIRST  DERIVATIVES  W.R.T. 

ARC  length  and  so  ON 

S1,S2,S3  DcNOTE  THE  IS T , 2N J , E T C , JE R IV  AT  I V ES  OF  THE 

ARC  length  H.R.T,  the  CENTRAL  ANGLL 
SO  CHANGE  OF  INTEGRATION  VARIABLE  FACTOR 

HTVLSI  calls  the  following  SUBPROGRAMS 

MVAL5  this  SJBROJTINE  USES  THE  QUANTITIES  IN  COMMON  BLOCK  4 TO 
COMPJTl  the  QUANTITIlS  IN  COMMON  BLOCKS  1 AND  2 
MDENSI  THIS  SUBROUTINE  USES  THE  VALUES  IN  COMMON  BLOCK  2 TO 

COMPJT:l  the  MATRIX  ENTRIES  FOR  OtTERMINING  THE  DENSITIES. 
THE  MATRIX  cNTRIES  ARE  STORED  IN  COMMON  BLOCK  5. 

SITER  this  subroutine  SOLVcS  THE  SYSTEM  OF  LINEAR 
EQUATIONS  FOR  THE  DENSITIES  BY  AN  ITERATION 
PROCESS,  this  ROUTINE  IS  APPLICABLE  TO  LINEAR 
SYSTEMS  OF  EQUATIONS  OF  A PARTICULAR  TYPE  ONLY. 

MODENS  THIS  SUBROUTINE  CONPUTES  THE  DERIVATIVES  OF  THE  DENSITIES 
MTVS  this  subroutine  COMPUTES  THE  CONTRIBUTIONS  TAA.TBB  AND  TCC 

TO  THE  TANGENTIAL  VELOCITIES. 

The  input  data  to  run  this  program  is  as  follows 

NOTE-THE  LETTER  8 DENOTES  A BLANK  IN  THE  INPUT  DATA. 

BBBNC 

BBBN(l) 

ETC. 

308N(NCI 


LLL  = 1 

PI=i.l4159?6&35!}979 
READ  900, NC 
NT  = 0 

DO  10  1=1, NC 
READ  900, N(I) 

10  NT=NT*N(I) 

REA J (5, NAM4) 

DO  20  U=1,NC 
REAO(5,900  M(U) 

M1P  = M(U)H'1 
M1M=M(U»-1 

READ(5,901)XX( U) ,YY(U) 
READ(5,U01)(A(I,U),I=1,M1P) 
20  R£AD(!>,  901)  (B(I,U»,I=1,M1M) 

call  MVALS 
CALL  MDENSI 
RlA J (5, NAM5) 

call  SITER(NT,K,EM,AAA,AA,L) 
IF  (L)  21,21,22 

21  PRINT  902, EM 
GO  TO  100 

22  WRITE (5  ,900) NT , K 
WRITE (5  ,901) EM 

WRITE (5  ,90 1)  (AA(I),I=1,nT) 
CALL  SITER(NT,K,EM,AAA,B0,L) 
IF(l)  2 3,23,  24 


23  PRIMT  902, 

G3  TO  100 

24  WRir£(5,900)NT,< 

WRIT-:  (5  ,901)  Ei 
WRIT:(5,901»  (3B( I) , I=1,NT) 

OJ  iO  J=1,NC 

CftLL  SITER(NT,<,aM,4AA,CC(l,J»,L> 

IF(L)  2B,26,27 
2b  PRI.TT  9 02, E'l 
GO  TO  100 

27  WRIT£(5  ,900) N( J)  ,< 

WRITE (b ,90 1) EM 

WRIT:<5,901) (3C(I,J) ,1=1, NT) 

00  CONTINUE 

WRITE(5  ,NAM5) 

call  mdoens 

WRITE (5  ,900) NT 

WRITE(5,901) (JAA(I),I=1,NT) 

writ: (5 ,900) NT 

WRITE  (5  ,901)  030(1)  , 1 = 1,  NT) 

nn  ^0  j=i,NC 

(5  ,90  0)  N(  J) 

40  E(5  ,901)  (DCC (I, J) ,I  = 1,NT) 

5  ,NAM7) 
w 1 T \IS 

rtRITE(  5 ,900) NT 

WRITE (5  ,90 1)  (T AA (I)  , 1 = 1, NT) 

WRIT£(5  ,900) NT 
WRITE(5,901) (T30(I) , 1=1, NT) 

DO  50  J=1,NC 
WRITE (5  ,900) N( J) 

50  WRITE  (5  ,901)  (TCCd,  J)  ,I=1,NT) 

WRITE (5  ,NAM3 ) 

100  STOP 

900  E0RMAT(-»I5) 

901  FOR1AT  ( 5E23.  !«♦) 

902  F3RMAT(*0EM=*£23.14,  •,  ITciRATION  PROCESS  DIO  NOT  CONVERGE*) 
END 


su'H3jt  in:  mocnsi 

COmN/  iT2/T  I(  2)  ,SO(  30)  ,S1  ( 30)  i S2(90)  ,S3  (80)  ,X  (dO)  , 

3 XK8(J)  ,X2(30),X3(30),Y(30),Y1(60),Y2(60),Y3(30),X1P(30), 
*♦  X2P(3C),X3P(30),Y1P(80),Y2P(80),Y3P(60) 

6 /3T4/C  (2) ,XX( 2) , YY { 2) ,M(2) ,A (51 ,2) ,8(49, 2) 

7 7JT5/AAA(30,80) ,IP(30) 

9 /-3T7/N  (2)  , NC,  NT 


c This  subroutine  computes  the  matrix  used  for 
c determining  The  source  densities. 


PI=3. 14153265353979 


100 


I?=ST-1 
00  iO  1=1,12 

ji=r«-i 

00  20  J=Jl,'iT 

U=(X(I)-X(J))*Y1P(U-(Y(I)-Y(J))*X1P(I) 
V=(X(II-X(J>)**2>(Y(I)-Y(J))**2 
20  Alft(  I,  J)  = ('J/V)  *S0(J) 

30  continue 

DO  50  1=2, NT 
J2=I-1 

00  •♦0  J = 1,J2 

U={X(I»-X(jn*YlP(I»-(Y(I)-Y(J))*XlP(I) 
V=(X(I)-X(J))**2*(Y(I)-Y(J))**2 
40  ftftAd,  J»  = (U/V»  *SO(J) 

50  continue 

00  oO  1=1, NT 

U=Y2(I)*X1P(I>-X2(I)*Y1P(I) 

V=X1(I) **2+Yl(I>**2 
60  AAA(I,I)=0. 5»(U/V)*30(I) 

RETURN 

ENO 


SUBROUTINE  SITER(N,K,£N, A,0,L) 
OINENSIJN  A(80,90),3(30),V(80),BC(dO) 


C THIS  SUBROUTINE  SOLVES  THE  LINEAR  SYSTEM  OF  EQUATIONS 

C (PI*I<-A)X=B  BY  AN  ITERATION  PROCESS.  HERE  PI  IS  AN  EIGENVALUE 

C OF  THE  MATRIX  A.  FIRST,  REWRITE  ThE  SYSTEM  OF  EQUATIONS  IN  THE 
C FORM  ( I^A/PI) X=B/PI  ANO  SET  AC=A/PI  ANO  B0=0/PI.  THE  MATRIX  AC 
C HAS  1 FOR  AN  EIGENVALUE.  NEXT  WE  WANT  TO  RtMOVE  FROM  THE  VECTOR 
C BO  any  eigenvectors  corresponding  to  the  eigenvalue  1.  THIS  IS 
C 30NE  BY  WRITING  X = ( 30/ 2 ) +-XC.  ONE  OBTAINS  THE  EQUATION  ( I 4-AC)  XC  =BC , 

C WHERE  3C= (B0-AC*30) /2, WHICH  MUST  BE  SOLVED  FOR  XC.  ONE  HAS  THEN 

C XC  = 3C-AC*  BC«-AC*(  AC*33) -♦ETC. , ANO  FINALLY  ONE  OBTAINS  X. 

C SITER  calls  the  FOLLOWING  SUBPROGRAM! 

C SCV  this  subroutine  GIVES  The  PRODUCT  OF  A 
C VECTOR  ANJ  A SCALAR 

C MATV  THIS  SUBROUTINE  GIVES  THE  PRODUCT  OF  A VECTOR 
C BY  A MATRIX  AND  DIVIDES  THE  RESULTING  VECTOR 

C BY  PI. 


PI  = J. 1415926535898 
C=1/PI 

call  SCV(N,C,3) 

DO  5 1=1, N 
5 V(I)=3(I) 

call  MATV(N,A,V) 

DO  10  1=1, N 
V(I»  =0. 5*(B(  I»  -V(in 
10  SC(I)=V(I) 

K=1 
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SI=-1 

15  call  MATV(N,A,V) 
iO  I=1,N 

^0  9:m  =g::(i)  vsi^vd) 

SI=-SI 
£*1=  3 

03  50  1=1, N 
U=A  5S(V  (II  ) 

►.  = A nxi  (Ei  , J) 

50  E5=; 

IE  (1.l-05-:lM)  **0,40,73 

4U  IP(?J-KI  -50,60, JO 
50  <=<♦! 

GO  TO  IJ 
60  L = 0 

GO  TO  90 
70  L = 1 

00  30  1=1, N 

30  8(1  I =0. 5»9 (I  H-3C( I» 

90  RETURN 
END 


REAL  FUNCTION  MFOU ( X , N , A , B , L ) 

Ol'IENSION  A(51,2»  ,3(49,2) 

»»«»  »*»«  ****  ****  **** 

C»,,*  *,*,  ,,,,  **,,  ,,,, 

C TRI3  PROGRAM  EVALUATES  THE  TRIGONOMETRIC  INTERPOLATING 

c polynomial  at  a given  point  X. 

C4***  *«4»  »♦*»  *♦»*  *■*■»»  ****  I*'4’** 

MFOU  = (A  (1,L)  4-A  (N+1,L  ) *COS  (N*X)  ) ♦O.S 
NM=N-1 
F = 0 

DO  15  K=1,NM 

15  F=F+A(K>1,L) •C0S(K*X)v3(K,L)*SIN(K*X> 

MF0J  = MF  OU4-F 

RETURN 

END 


real  function  MOFOJ(X,N, A,8,L) 
dimension  A(51 ,2) ,3(49,2) 

C,,,,  ♦*,,  ,,,,  ,,»*  ,,,,  *,,,  ,,,,  ,,,, 

Q*»»*  »*»♦  ♦«,»♦  »♦♦♦  •♦♦♦  ****  ♦♦♦» 

C THIS  program  EVALUATES  THE  DERIVATIVE  OF  THE  INTERPOLATING 

C POLYNOMIAL  AT  A GIVEN  VALUE  X. 

Q,,,,  ,,,♦  ,,,,  ,,,,  ,,,,  ,,,,  ,,,,  ,4,*  44,4 

NM=N-1 

NP=N«-1 

F=-.5*N*A(NP,L) *SIN(N*X) 
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00  15  K=1,N“1 

15  F = r-<*A(K4-l,L)»SIN(K*X)^K*3(K,U*C0S(K*X) 

MJFO J=F 
return 
EMD 

REAL  FUMCTION  H02F0U ( X , N , A , 3, L) 

DIMENSION  A(51 , 2)  ,3(  (♦9,2) 


This  proG'JAm  -»/aluates  the  second 

OERIy/ATIy/E  OF  THE  INTERPOLATING 
POLYNOMIAL  AT  A GH/EN  y/ALUE  X, 


NM=N-1 
NP  = t01 

F=-.5*N*N*A(N3,L»*C0S<N*X» 

00  15  K=1,NM 

IK  F=F-<*K*(A (<tl,L)*COS(K*X) ♦9(K,L)*3IN(K*X) ) 
M02F0U=F 
RETJRN 
ENO 


REAL  FUNCTION  MO 3F0U ( X , N , A , B, L ) 
DI1ENSI0N  A(51, 2) ,3(49,2) 


C this  program  evaluates  the  THIRD  OERI- 
C VATIVt  OF  THE  INTERPOLATING  POLYNOMIAL 
C AT  A GIVEN  VALUE  OF  X 


NM=N-1 

NP=N*-1 

F=.5*(N**i)*A{NP,L)*SIN(N*X) 

00  15  K=1,NM 

15  F = Ff(K**3)*(A(<4-l,L)*SIN(K*X)-9(K,L)  *C0S(t<*X)  ) 
M03F0U=F 
RETURN 
END 


103 


PROGRAM  KATRCFAPES, INPUT, OUTPUT) 

CO'i'ION/Kl/XA  (60,2),Xa(50,2)/K2/TH(70),RO(73),R2(70) 
NAH£LIST/NAH/LLL/NAN  1/LlL/NA'12/LLU/NAM3/LLL/NAMh/LLL 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


THIS  15  th;  control  program  for  transforming  an 

AIRFOIL  PROFILE  IN  TO  A NEAR  CIRCLE  RY  MEANS  OF 
A KARMAN  TREFFETZ  LIKE  TRANSFORMATION. 

THE  ENJ  PRODUCTS  OF  THIS  PROGRAM  ARE  THE 
COOROI TATES,  IN  POLAR  FORM  OF  THE  IMAGES 
OF  THE  profile  DATA  ON  THE  NEAR  CIRCLE 
AND  THE  SECOND  OERI>/ATIVES  AT  THESE  POINTS 
OF  The  function  defined  9y  These  near 

CIRCLE  POINTS.  THESE  QUANTITIES  DETERMINE 
THE  CUaiC  SPLINE  FIT  OF  THE  NEAR  CIRCLE. 

the  result  from  This  program  are  written  on 
A file  which  the  user  may  name  and  store. 

NA(NB)  - the  number  OF  DATA  POINTS  ON  THE  UPPER (LOWER) SURFACE. 
XA(Xa>  - THE  SET  OF  COORDINATES  OF  THE  DATA  POINTS  ON  THE 

upper(lower)  surface  of  the  airfoil,  the  coordinates 

OF  the  LEADING  AND  TRAILING  EDGES  ARE 
INCLUDED  IN  DOTH  SlTS.  THESE  SETS  ARE  INDEXED 
beginning  at  the  leading  EDGL. 

TH(I)  - the  argument  OF  THE  NEAR  CIRCLE  DATA 
RO(I)  - THE  MAGNITUDE  OF  ThE  NEAR  CIRCLE  DATA 
R2(I)  - THE  SECOND  DERIVATIVE  OF  THE  CUBIC 

SPLINE  FIT  OF  the  NEAR  CIRCLE  DATA. 

TA(TB)  - THE  ANG^E  THE  UPPER(LOWER)  SURFACE  MAKES  THE  X-AXIS 
AT  THE  TRAILING  c-DGE. 

G - THE  EXPONENT  USED  IN  THE  SUB- 
PROGRAM GP. 

T - THE  IMAGE  IN  THE  WI-PLANE  OF  THE  POINT 
AT  infinity  in  THE  Z-PLANE  IS  EXP(IT). 

THIS  PROGRAM  CALLS  THE  FOLLOWING  SUBPROGRAMS 
FlI  computes  W=  ( 14-Z)  / (1-Z)  . W IS  IN  POLAR  FORM 
AND  Z IS  A DATA  POINT  ON  THE  SURFACE. 

REC  CHAN^^ES  A SET  OF  POINTS  IN  POLAR  FORM  TO 
rectangular  form. 

GP  COMPUTES  Wl= (W*EXP(-I*TA) ) **G.  W1  IS  IN 
POLA^FORM.  G=pI/(TB-TA) . 

FL2  COMPJTES  Wj  - ( W 1 ♦ T ) / ( W 1- T ) . W3  IS  IN  POLAR  FORM  AND  T 
IS  The  iMAGt.  IN  THE  wi  plane  of  The  point  at  infinity 
IN  THE  Z-PuANE. 

POLAR  CONVERTS  R_C T ANGULAR  COORDI4ATES  TO  POLAR  COORDINATES. 
ICOF  GIV-^S  A CUBIC  SPlINE  FIT  OF  THi.  RADIUS  OF  THE  NEAR 
circle  as  a funtion  of  The  polar  angle. 

ICOF  calls  LUJ. 

LU3  A routine  FOR  HE  DIRECT  SOLUTION  OF  SYSTEMS 
DF  LINEAR  ^QUATIONj  TAKEN  FROM  THE  ARL 
LINEAR  algebra  uIBRARY.  IT  INCLUDES  LEQS3. 

THE  NiAR  CIRCL:l  IMAGES  OF  THE  PROFILE  DATA  POINTS  ARE  INDEXED 
FROM  1 TO  K=NA*NB-1  WITH  RO(K)=RO(l)  AND  TH ( K) = TH ( 1 )♦ 2PI , 


Li.L=  1 

PI  = J.  14  1592oE«ib398 
REAJ  90J,Na,Nn 

RtA3  901,((Xft(I,J),J=l,2),I=l,NA) 
REAJ  901,((XB(I,J),J=1,2),I=1,NB) 
WRITE (5 ,900) Nft , NB 

WRITE(5,g02)((X4(I,J),J=l,2),I=l,Nft) 
WRITE  (5  ,902)  nX9(I,J),J=l,2),I  = l,Na) 
WRITE (5  ,NAM) 

WRITE(5  ,NAM1) 

UA  = XA (N  a-1 , 1 ) ♦ 1 . 

VA=XA(NA-1,2) 

U9=X3(N9-l,l)f 1. 

V3  = X9(N  )-l,2) 

CALL  POLAR  (UA,  y/A,RA  , TA) 

CALL  POLAR (UB, VB,R9, TB) 

G = PI/ ( T 3-T A) 

call  FLl(iJA,N9) 

C REC  IS  CALLED  TO  MORMALIZE  ThE  RESULTS 
C FROM  PLl. 

call  R£C(MA,N3) 
call  ^,P  (NA,MB,  j,TA) 

CALL  REC(MA,N3) 

T=G* (PI -TA) 
call  FL2(NA,NB,T) 

K = “JA<-NB-1 
TH(1) =0  . 

TH(<) =2  .’PI 
CALi-  ICOF(K) 

WRITE (5 ,90  3) TA , TB,G, T 
WRITE (5  ,NAM2 ) 

WRITE (5  ,900) K 
WRITE(5  ,903)  ( T H ( I ) , I = 1 , K ) 

WRITE (5  ,90  3)  <RO ( I ) , I =1 , K ) 
WRITE(5,903) (R2(I» ,I=1,K) 

WRITE (5  ,NAM3) 

STOP 

FORMAT ( hIS) 

FORMATtaElO.O) 

format  < bEI’. 8) 

FORMAT ( 9t23. 14) 

ENO 


900 

9C1 

902 

903 


SJBRO'JT  INE  FLI  (N,M) 

common/ <1/ XA (p  0 ,2) , X 3(50 , 2) /K2/TH( 70 ) ,R0( 70)  ,R2  (70 ) 


C THIS  PROGRAM  COMPUTES  W*(l»2)/H-Z)  AT  A GIVEN  SET 
C OF  POINTS  Z.  Z IS  COMPLEX. 

C W IS  IN  polar  FORM, THAT  IS, 

C M(J)=RO( J) •EXP(I*TH(J) » . W IS  ^Ht 
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C L-iA&E  JF  THE  POINTS  ON  ThE  'JPPER  AND 
C LOWER  PROFIL:.  SURFACES.  THESE  IMAGE 
C POINTS  ARE  INJlX:.0  CONS  tCUT  I VEL  Y 
C AROUND  the  PROFILE  CONTOUR,  RESINNING 
C AT  THE  LEADING  EDGE  OF  THE  UPPER  SURFACE. 
C FLI  calls  THE  SUBPROGRAM  POLAR. 


NM=H-1 

00  10  1=1, NM 
UA=-XA ( I, 1) -1 
VA  = -XA ( 1,2 ) 

U3=-UA- 2 
y/B=-VA 

call  P0LAR(UAWA,K4,TA) 
call  P0lAR(UB,\/J,R3,T3) 
R0(I)  =RA/RS 
10  TH(I)=TA-T3 
RO(N) =0 
Th(N) =C 
MM=M-1 
NP  = N+  1 
NPH=N-t-HM 
J = 0 

00  20  I=NP,NPM 
J = J*1 

UA=-X3(H-J,1)-1 

VA=-X3(M-J,2) 

Ua=-UA- 2 
\/B  = -ifA 

call  P0LAR(UA, VA,RA, TA) 
call  POLARTUB, V3,RB, TO) 
RO ( I)  =RA/RR 
20  TH(I)=TA-TB 
RETURN 
END 


SJIROUTINE  GP(N,M,G,TA) 

COM  ION/ <2/ TH (7  0 > ,RO(70) ,R2(70) 


C COMPUTES  W 1=( W*EXP(-I*TA» » »*G.  W1  IS  IN  RECTANGULAR  FORM. 

c GP  calls  polar. 


J=1 

NM=N-1 

5 00  10  I=J,NM 

call  POlAR<RO(I) ,TH(I) ,R,T) 
R0(I) =R*»G 
10  TH(I)  =(  r-TA) 
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IF(J-l)  15,15,20 
15  JsN*! 

NM  = N»M-  1 
GO  TO  5 
20  RETU'^N 

END 


SJDROJTINE  'tEC(N,M» 

COMMON/ K2/TH (70 ) , RO ( 7 0 » , R2 ( 7 0 ) 


C this  program  converts  a set  of  points  in  POLAR 
C FORM  TO  A SET  OF  POINTS  KITH  RECTANGULAR 
C COORDINATES.  THE  NEW  COORDINATES  ARE  STORtO 
C ON  TOP  OF  THE  OLD. 


K=N+M-1 
DO  10  1=1, K 
R=RJ(I) 

B=COS (T  H(I ) ) 
C = SIN(TH(I)  ) 
RO(I) =R*B 
10  TH(I>=R*C 
RETURN 
END 


SUBROUTINE  FL2(N,M,T) 

COMMON/ <2/ TH (70 ) , RO ( 70 ) , R2 ( 70 ) 


C THIS  PROGRAM  COMPUTES  W 5=  (HK-T  ) / ( Wl-T  ) . W3  IS  IN 
C POLAR  FORM  AND  T IS  THE  IMAGE  IN  ThE  Wl-PLANE  OF 
C THE  POINT  AT  INFINITY  IN  THE  Z-PLANE. 


K=N*M-1 
U=COS  (T  ) 

V = SIN(T  ) 

DO  10  1=1, K 
UA  = RO  (I  H-U 
VA=TH(I )+V 
U3=R0(H-J 
V0=TH(I)-V 

CALL  POLAR(UA, VA,RA, TA) 
CALL  POLAR(UB, V3,RU, TH) 
R=RA/R9 
TT  = TA-T  J 
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X = R*GOS  (TT  » 
Y=R*SIN(TT) 

C4Li-  POLA^<x, ir, RO<n 

(RETURN 

ENU 


SJ-3'iOUTINE  POLftR(X,Y,R,T) 


C THIS  CONVERTS  RECTANGULAR  COORJINATES  TO 

C POLAR  COORJINATE5, 


• ♦ • 


5 

1 

8 

9 

10 
11 
15 


20 

21 

22 

25 

26 

27 

3u 

31 

4C 


S=3. 1415925535898 
A=X*  Y 

IF(A)  15,5,15 
IF(X»  10,5,11 
IF(Y)  8,7,9 
R = 0 
T = 0 

GO  TO  40 
R=A9S(Y) 

T=(3.*S)/2. 

GO  TO  40 
R = Y 

T=S/2. 

GO  TO  40 
R = A3S (X  > 

T = S 

GO  TO  40 
R=X 
T = 3. 

GO  TO  40 
R=SaRT ( X*XfY*Y  > 

U=A9S(Y/X» 

V = AT  AN ( U) 

IF  (d)  20,30,25 

IF(X)  21,30,22 
T = S-V 
GO  TO  LO 
T=2.*5- V 
GO  TO  40 
IF  (X)  26,30,27 

T = S*V 
GO  TO  40 
T = V 

GO  TO  40 
PRINT  31 

FOR1AT(*05U9R0 JTINE  POLAR  PROUUCT  XY  IS  BOTH  ZERO  AND 
INOT  ZERO  »> 

STOP 

RETURN 
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ENO 


SUHOUTINE  IC3F(N) 

OIMINSION  ft(70,70),IP(70) 

CO -1-1  ON/  <^/  X(7u)  ,Y(7!)),Z(70) 

NA'itLIS  T/NA'I/LLL/NA'i  l/LLL/NAi2/LLL/NAM3/LLL/NAM4/LLL 


c this  p-^ogram  determines  a cjbic  spline  fit  for  a periodic 

C FUNCTION.  the  FORMULA  FOR  DETERMINING  ThE  VALUE  OF  THE 
C FUNCTION  AT  A POINT  T DEPENDS  UPON  THE  FUNCTION  VALUES  AT 
C the  ENDPOINTS  OF  THE  INTERVAL  ON  WHICH  T LIES  AND  THE 
C VALUES  OF  the  2ND  DERIVATIVES  OF  THE  FUNCTION  AT  THESE 

c ENOPoius.  This  program  determines  the  values  of  the  2nd 

C DERIVATIVES  jF  ThE  FUNCTION  AT  THE  NODES  WHERE  THE 
C FUNCTION  values  ARE  GIVEN. 

C ICOF  CALLS  THE  SUBPROGRAM  LU3 
C LU3  IS  A ROJTINl  FOR  THE  DIRECT  SOLUTION  OF 
C OF  SYSTEMS  OF  LINEAR  EQUATIONS. 


LlL=1 

PI=  3. 14  1592oS35398 
NM=N-1 

DO  10  1=1, NM 
DD  6 J=1,NM 
a A(I,J)  = 0 

10  CONTINUE 

K = N-3 

OJ  15  1=1, < 

J = I 

A(I,J»=X(I41)-X(I) 

A (I , J41 )=2*( X( I+2)-X  (I)  ) 

A (I, Jf2  »=X (I *2 ) -X (I f 1) 

15  Z(I)=6*(<Y(If2)-Y(I<-in/A(I,If2»-(Y(l41)-Y(I)>/A<I,I)) 

A (N-2, 1 )=X (N) - X ( N-1 » 

A IN-2 ,N-2»  =X (N-l)-X (N-2) 

A( N-2, N-1)  = 2*(X(N)-X(N-2)) 

Z (N-2) =o* ( (Y (1 ) -Y (NM) ) /A (N-2, 1 )- ( Y(NM) -Y (N-2) ) / A(N-2,N-?) ) 
A(NM,1)=2*(A(N-?,1)>A(1,1)) 

A (NM, 2) =A ( 1, 1) 

A (NM, NM ) =A (N-2, 1) 

Z(NM)=b*((Y(2)-Y(l))/A(l,l)-(Y(l)-Y(NM))/A(N-2,D) 

CALL  LUi(NM,A,70,IP) 

CALc  LE aS3(NM,A,70,IP,Z) 

Z(N)=Z(  1) 

RETURN 

900  F0RMAT(4I5) 

903  FORMAT ( 5E23.  14) 

END 
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